東 海 大 学 紀 要 産 業 工 学 部 5( 年 )7 頁 3 頁 Bll. School o Idstrl Egeerg ok U. 5() pp.7-3 LDU 分 解 時 における 固 有 値 についての 行 列 式 の 微 分 柏 木 光 博 Derte o Determt wth respect to Egele LDU Decomposto b Mtshro KASHIWAGI (Receed: September 3 Accepted: Febrr 3) Abstrct hs pper presets derte o determt wth respect to egele LDU decomposto d the mercl lgorthm. he proposed method s eecte especll the bd mtr tht c se storge cpct shrpl compred wth the desed mtr. hs method comptes the sglr pot oler lss d egele problem sg proposed method. Otpts o mercl emples show tht the proposed de works well d ldt throgh mercl smlto. Ke Words : Derte o determt LDU decomposto Bd mtr No-smmetrc mtr *. はじめに 著 者 は 共 役 勾 配 法 の 理 論 の 中 の 係 数 行 列 の 逆 行 列 算 定 式 とトレース 理 論 を 使 って 固 有 値 問 題 の 行 列 式 の 固 有 値 による 微 分 式 を 提 案 した ). 特 に 固 有 値 解 析 法 に 応 用 し 効 果 的 であることを 示 してきた -6).しかしそこ での 提 案 式 は 共 役 勾 配 法 などの 反 復 法 用 であり LDU 分 解 などの 直 接 法 には 使 えない.ここでは 共 役 勾 配 法 に 関 する 式 ) を 参 考 にして 直 接 法 の 代 表 的 解 法 である LDU 分 解 時 において 非 対 称 行 列 の 固 有 値 についての 行 列 式 の 微 分 を 示 す.またニュートン ラフソン 法 にお ける 解 は 行 列 式 の 微 分 によって 変 化 するのでそれを 増 分 量 として 利 用 できる. 増 分 量 の 設 定 が 自 動 化 できるの で 非 線 形 解 析 への 影 響 が 大 きい 基 本 的 パラメーターを 減 らすことができる. 提 案 法 は 密 行 列 の LDU 分 解 や 密 行 列 に 比 べ 記 憶 容 量 を 大 幅 に 節 約 できる 帯 行 列 の LDU 分 解 に 対 してもバンド 幅 内 の 要 素 を 使 って 計 算 できるの で 密 行 列 の 計 算 をあたかも 帯 行 列 の 計 算 をするかの ようにできる. 提 案 式 を 含 まない 帯 行 列 の 計 算 時 間 と 比 べてもそれほどの 増 加 にならない. 以 下 にその 理 論 と 詳 細 なアルゴリズムおよび 移 流 拡 散 方 程 式 を 離 散 化 した * 産 業 工 学 部 教 授 非 対 称 行 列 による 数 値 実 験 例 を 示 す.. LDU 分 解 時 における 固 有 値 による 行 列 式 の 微 分 式 固 有 値 問 題 は 次 のように 表 される. A を の 非 対 称 行 列 とすると 標 準 固 有 値 問 題 は (.) として 示 される.ここに は 固 有 値 を は 固 有 ベク トルを 表 している. 式 (.)が 自 明 な 解 をもつためには が 特 異 でなければならないので det[ ] (.) ここで を 以 下 のようにおく. det[ (.3) ] トレース 理 論 によって ' trce trce (.4) - 7 -
LDU 分 解 時 における 固 有 値 についての 行 列 式 の 微 分 A sm. [ ] (.5) とおく. LDU 分 解 時 には A LDU (.6) になる.ここに L [ ] (.7) ( ) D [ d U [ L d d ] ] d ( ) L を 以 下 のようにおく. [ g LL I g ] g g ( ) (.8) (.9) (.) ( I は 単 位 行 列 )を 展 開 して 整 理 すると U U g g g (.) k k 次 に U を 以 下 のようにおく. [ h h h ] h ( ) U I を 展 開 して 整 理 すると h k (.) h h (.3) k k k 式 (.)と(.3)から g と h は 全 要 素 について 計 算 する 必 要 があるが 半 帯 幅 外 の 要 素 については お よび となるので 半 帯 幅 内 の と を 使 って 計 算 できる.また 半 帯 幅 が 狭 いほど 計 算 時 間 は 短 くなる. 行 列 式 の 微 分 は 式 (.4)より 式 (.6)の 逆 行 列 の 対 角 要 素 だけを 求 めて 加 算 すればよい. 式 (.7)(.8)(.9) (.)(.)から U D L 求 めると 式 (.4)は 以 下 のようになる. ' d k g k h d kk k を 展 開 し 対 角 要 素 の 和 を (.4) 以 上 から 行 列 式 の 微 分 には LDU 分 解 時 の D 行 列 と L 行 列 およびU 行 列 の 逆 行 列 の 要 素 から 成 っているこ とが 分 かる. 上 述 のようにバンド 幅 内 の 要 素 を 使 って 計 算 できるので 密 行 列 の 計 算 をあたかも 帯 行 列 の 計 算 をするかのようにできる.よって 提 案 式 を 含 まない 帯 行 列 の 計 算 時 間 と 比 べてもそれほどの 増 加 にならな いと 予 想 される. LDU 分 解 プログラムに 補 足 したプログラム(ベクト ルを 本 追 加 )により 簡 単 に ' が 得 られる.この 値 はニュートン ラフソン 法 やデュラン カーナー 法 な どによる 数 値 解 析 ( 固 有 値 に 近 づくと () は に 近 づく が 実 際 の 計 算 では ' を 使 うので 問 題 ない)に 便 利 に 利 用 できる.また 提 案 法 は 連 立 一 次 方 程 式 を 解 く 過 程 で 簡 単 に 求 めることができる 方 法 となっている. 3. 提 案 法 のアルゴリズム 3-. 密 行 列 用 アルゴリズム 最 初 に 密 行 列 に 対 するアルゴリズムを 示 す. 配 列 や 変 数 は 以 下 による ) LDU 分 解 と 行 列 式 の 固 有 値 による 微 分 値 の 算 定 入 力 データ A : ge coecet mtr -dmeso rr s A() b : work ector -dmeso rr s b() c : work ector -dmeso rr s c() : ge order o mtr A d ector b eps : prmeter to check sglrt o the mtr otpt 出 力 データ A : L mtr D mtr d U mtr - 8 -
柏 木 光 博 -dmeso rr s A() d : derte o determt err : error code = or orml eecto = or sglrt 3LDU 分 解 do = <d()> do k=- A()=A()-A()*A()*A() (bs(a())<eps) the err= retr ed <l() d ()> do =+ do k=- A()=A()-A(k)*A(kk)*A(k) A()=A()-A(k)*A(kk)*A(k) A()=A()/A() A()=A()/A() err= 4 行 列 式 の 固 有 値 による 微 分 d= <()>. do = d=d-/a() <()> do =- do =+ b()=-a() c()=-a() do k=-- b()=b()-a(+k)*b(+k) c()=c()-a(+k)*c(+k) d=d-b()*c()/a() ) 解 の 算 定 入 力 データ A : L mtr D mtr d U mtr -dmeso rr s A() b : ge rght hd sde ector -dmeso rr s b() : ge order o mtr A d ector b 出 力 データ b : work d solto ector -dmeso rr 3 前 進 代 入 do = do =+ b()=b()-a()*b() 4 後 退 代 入 do = b()=b()/a() do = =-+ do =- b()=b()-a()*b() 3-.バンド 行 列 用 アルゴリズム 次 にバンド 行 列 に 対 するアルゴリズムを 示 す. 以 下 に 帯 行 列 の 概 要 を 図 - に 示 す.ここに は 左 バンド 幅 であり は 右 バンド 幅 である.それぞれに 対 角 要 素 は 含 まないので 帯 幅 は 図 - 帯 行 列 と 帯 幅 である. ) LDU 分 解 と 行 列 式 の 固 有 値 による 微 分 値 の 算 定 入 力 データ A : ge coecet bd mtr -dmeso rr s A(l++) b : work ector -dmeso rr s b() c : work ector -dmeso rr s c() : ge order o mtr A l : ge let hl bd wdth o mtr A : ge rght hl bd wdth o mtr A eps : prmeter to check sglrt o the mtr - 9 -
LDU 分 解 時 における 固 有 値 についての 行 列 式 の 微 分 出 力 データ A : L mtr D mtr d U mtr -dmeso rr s A(l++) d : derte o determt err : error code = or orml eecto = or sglrt 3LDU 分 解 do = <d()> do =m(-m(l))- A(l+)=A(l+)-A(l+-(-))* A(l+)*A(l++(-)) (bs(a(l+))<eps) the err= retr ed <l()> do =+m(+l) sl=a(l+-(-)) do k=m(-l)- sl=sl- A(l+-(-k))*A(kl+) *A(kl++(-k)) A(l+-(-))=sl/A(l+) <()> do =+m(+) s=a(l++(-)) do k=m(-)- s=s-a(l+-(-k))*a(kl+) *A(kl++(-k)) A(l++(-))=s/A(l+) err= 4 行 列 式 の 固 有 値 による 微 分 d= <()> do = d=d-/a(q) <()> do =- do =+m(+l) b()=-(l+-(-)) do k=-- b()=b()-(l+-(-)+k)*b(+k) do =+l+ b()=.d do k=l b()=b()-(k)*b(-l-+k) do =+m(+) c()=-(l++(-)) do k=-- c()=c()-(+kl++(-)-k)*c(+k) do =++ c()=.d do k= c()=c()-(l++k)*c(-l-+k) do =+ d=d-b()*c()/(l+) ) 解 の 算 定 入 力 データ A : ge decomposed coecet bd mtr -dmeso rr s A(l++) b : ge rght hd sde ector -dmeso rr s b() : ge order o mtr A d ector b l : ge let hl bd wdth o mtr A : ge rght hl bd wdth o mtr A 出 力 データ b : solto ector -dmeso rr 3 前 進 代 入 do = do =m(-l)- b()=b()-a(+-(-))*b() 4 後 退 代 入 do = =-+ b()=b()/(l+) do =+m(+) b()=b()-(l++(-))*b() - -
柏 木 光 博 4. 数 値 実 験 提 案 式 の 有 効 性 を 示 すために 非 対 称 行 列 を 有 する 標 準 固 有 値 問 題 に 関 する 例 題 による 数 値 実 験 を 行 った.こ こでの 数 値 実 験 はすべて 倍 精 度 演 算 とし 求 める 固 有 解 の 個 数 は 個 とした.アルゴリズムと 数 値 実 験 結 果 を 以 下 に 示 す. 計 算 には 今 回 はスーパーコンピューター を 対 象 としていないことおよび 最 近 のパーソナルコン ピューターの 性 能 向 上 などを 考 慮 し Itel(R) Core 7 3.GHz RAM. GB Wdows 7 PGI Fortr Workstto (pgws64-6) を 使 用 した. 以 下 に 詳 細 を 記 す. 4-.アルゴリズム ここでは 逆 ベキ 乗 法 とシフト 逆 ベキ 乗 法 の 提 案 式 用 ア ルゴリズムの 詳 細 を 示 す.これらのアルゴリズムにより 数 値 実 験 を 行 った. A は 行 列 は 固 有 値 は 対 応 する 右 固 有 ベ クトル は 対 応 する 左 固 有 ベクトル を 原 点 移 動 量 とすると 標 準 固 有 値 問 題 はA I およびA I のように 表 せる. 4--. 逆 ベキ 乗 法 () として A I ゴリズムでは 常 に である. () 初 期 ベクトル () と () の LDU 分 解 を 行 う. 本 アル を 一 様 乱 数 によって 生 成 する. 最 大 反 復 回 数 を 設 定 する.グラム シュミットの 双 直 交 化 法 により 既 に 得 られた 個 の 右 固 有 ベクトル と 左 固 有 ベクトル の 抜 き 取 りと 正 規 化 および 初 期 近 似 固 有 値 を 次 式 で 求 める. () () () ( / ) () () () () ( ) () () () ( / ) () () () () () ( ) () 逆 ベキ 乗 法 による 反 復 計 算 を 行 う. 反 復 時 での () と () ( ) ( ) と A I () A () ( ) ( ) は A I によって 求 める.ここに は 転 置 行 列 の 逆 行 列 を 表 している.グラム シ ュミットの 双 直 交 化 法 により 既 に 得 られた 固 有 ベ クトルの 抜 き 取 りと 正 規 化 および 近 似 固 有 値 を 次 式 で 求 める. ( ) ( ) ( ) ( / ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( / ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) A 収 束 判 定 を 行 う. 最 大 反 復 回 数 ( 本 数 値 実 験 では m )を 超 えたら 計 算 をストップ. と を 収 束 判 定 値 とすると ( ) ( ) ( ) ( 本 数 値 実 験 では 4 ) あるいは 残 差 ベクトル: r 計 算 し r A ( ) ( ) A ( ) ( ) ( ) ( ) ( 本 数 値 実 験 では 8 ) であれば 収 束. 次 の 固 有 値 を 求 める 場 合 ()に 戻 るが 求 めない 場 合 (3)へ. 収 束 していない 場 合 ()のに 戻 る. (3) 最 後 に 求 められた 全 ての 固 有 値 と 固 有 ベクトルに 対 して 固 有 値 と 固 定 原 点 移 動 量 との 差 の 絶 対 値 が 小 さい 方 から 順 に 並 び 替 えを 行 う. 4--. シフト 逆 ベキ 乗 法 () として A I の LDU 分 解 を 行 う. () 初 期 ベクトル () と () を を 一 様 乱 数 によって 生 成 する. 最 大 反 復 回 数 を 設 定 する.グラム シュミットの 双 直 交 化 法 により 既 に 得 られた 個 の 右 固 有 ベクトル と 左 固 有 ベクトル の 抜 き 取 りと 正 規 化 および 初 期 近 似 固 有 値 を 次 式 で 求 める. () () () ( / ) () () () () ( ) () () () ( / ) () () () () () () () ( ) A () シフト 逆 ベキ 乗 法 による 反 復 計 算 を 行 う. 反 復 時 での () と () ( ) ( ) と A I ( ) ( ) は A I によって 求 める.ここに は 転 置 行 列 の 逆 行 列 を 表 している.グラム シ ュミットの 双 直 交 化 法 により 既 に 得 られた 固 有 ベクトルの 抜 き 取 りと 正 規 化 および 近 似 固 有 値 を 次 式 で 求 める. - -
CPU tme (sec) LDU 分 解 時 における 固 有 値 についての 行 列 式 の 微 分 ( ) ( ) ( ) ( ) ( / ( ) ( ) ( ) ( ) ( ) ( ) ) ( ) ( ) ( / ) ( ) ( ) ( ) ) ( ) ( ( ) ( ) ある 反 復 回 数 毎 (ここでは 反 復 回 数 )に 原 点 ( ) 移 動 量 を と 変 更 し I の LDU 分 解 A を 行 う.の 計 算 を 実 行 する. 次 に 今 回 求 めた FD( ' ) と 前 回 求 めた FD が 同 符 号 なら 原 点 移 ( ) 動 量 を / FD と 変 更 し I の A LDU 分 解 を 行 う.の 計 算 を 実 行 する. A 3 収 束 判 定 を 行 う. 最 大 反 復 回 数 ( 本 数 値 実 験 では m )を 超 えたら 計 算 をストップ. と を 収 束 判 定 値 とすると ( ) ( ) ( ) ( 本 数 値 実 験 では 4 ) あるいは 残 差 ベクトル: r 計 算 し r A ( ) ( ) A ( ) ( ) ( ) ( ) ( 本 数 値 実 験 では 8 ) を l l l l l l l l l l l m 上 記 以 外 のl m l Ω P( ) Ω であれば 収 束. 次 の 固 有 値 を 求 める 場 合 ()に 戻 るが 求 めない 場 合 (3)へ. 収 束 していない 場 合 ()のに 戻 る. (3) 最 後 に 求 められた 全 ての 固 有 値 と 固 有 ベクトルに 対 して 固 有 値 と 固 定 原 点 移 動 量 との 差 の 絶 対 値 が 小 さい 方 から 順 に 並 び 替 えを 行 う. 4-. 標 準 固 有 値 問 題 ( 移 流 拡 散 方 程 式 ) ここでは 定 常 問 題 における 移 流 拡 散 方 程 式 を 離 散 化 して 生 じる 係 数 マトリックスを 扱 った. 多 くの 極 近 接 根 をもつ 求 めにくい 対 象 である. 自 由 度 と 原 点 移 動 量 を 変 化 させながら 数 値 実 験 した. 以 下 のような 次 元 モデル を 考 える( 図 -). 境 界 条 件 は 境 界 値 を とする. 移 流 拡 散 方 程 式 に 中 心 差 分 を 用 いて 離 散 化 を 行 うと A の 固 有 値 問 題 が 得 られる.ここで 係 数 行 列 A の 各 要 素 は 次 式 となる. l A k l とおけば 図 - 次 元 モデルの 一 例 表 - 次 元 モデルのパラメーター tpe l l A 5 5 5 5 5 B C 5 5 5 5 5 D 4 8 Ierse power method 6 Sht erse power method 4 5 5 4 Nmbers o reedom 図 -3 次 元 モデルの 計 算 時 間 l - -
柏 木 光 博 ここに はそれぞれ 方 向 の 分 割 幅 である. ここでは とした.セルペ クレ 数 は で 行 列 A は 弱 対 角 優 位 行 列 となる. 数 値 実 験 に 用 いたパラメータ を 表 - [5] 柏 木 光 博 一 般 固 有 値 問 題 における 前 処 理 付 共 役 勾 配 法 による 大 次 元 スパース 対 称 行 列 の 固 有 解 の 一 算 定 法 日 本 建 築 学 会 構 造 工 学 論 文 集 Vol.56B445-45.3 [6] 柏 木 光 博 一 般 固 有 値 問 題 におけるランチョス 逆 ベキ 乗 法 による スパース 対 称 行 列 の 中 間 固 有 対 の 一 算 定 法 日 本 建 築 学 会 構 造 系 論 文 集 Vol.6648-88.6 に 示 す. 図 -3 は 表 - のモデルAの4タイプについて 逆 ベキ 乗 法 とシフト 逆 ベキ 乗 法 の 計 算 時 間 を 描 いている. 逆 ベ キ 乗 法 ではタイプ A と B については 解 析 できたがタ イプ C と D については 解 析 できず 次 元 数 が 増 加 すると 解 析 不 能 となった.しかし 提 案 のアルゴリズムによる シフト 逆 ベキ 乗 法 では4タイプ 共 解 析 が 可 能 であった. 提 案 法 は 計 算 時 間 や 反 復 回 数 も 少 なく 安 定 した 解 法 とな っている. 5. 結 び 直 接 法 の 代 表 的 解 法 である LDU 分 解 時 に 同 時 に 非 対 称 行 列 の 固 有 値 による 行 列 式 の 微 分 を 求 められることを 示 した.またニュートン ラフソン 法 における 解 は 行 列 式 の 微 分 によって 変 化 するのでそれを 増 分 量 として 利 用 することが 可 能 である. 増 分 量 の 設 定 が 自 動 化 でき るので 非 線 形 解 析 への 影 響 が 大 きい 基 本 的 パラメータ ーを 減 らすことができる. 提 案 法 は LDU 分 解 を 対 象 と しており 帯 行 列 は 密 行 列 の LDU 分 解 に 比 べ 記 憶 容 量 を 大 幅 に 節 約 できる.LDU 分 解 プログラムに 補 足 した プログラム(ベクトルを 本 追 加 )により 簡 単 に 固 有 値 による 行 列 式 の 微 分 が 得 られる.この 値 はニュート ン ラフソン 法 やデュラン カーナー 法 などによる 数 値 解 析 に 便 利 に 利 用 できる.また 提 案 法 は 連 立 一 次 方 程 式 を 解 く 過 程 で 簡 単 に 求 めることができる 方 法 であり 非 線 形 問 題 や 固 有 値 問 題 などに 対 して 特 に 有 効 な 方 法 と 考 えられる. 参 考 文 献 [] 柏 木 光 博 共 役 勾 配 法 による 最 大 あるいは 最 小 固 有 解 の 一 算 定 法 日 本 計 算 工 学 会 論 文 集 Vol.-5999.5 [] 柏 木 光 博 共 役 勾 配 法 による 大 次 元 スパース 対 称 行 列 の 固 有 解 日 本 応 用 数 理 学 会 論 文 誌 Vol.5 9-435.3 [3] KshwgM. A Method or Determg Egesoltos o Lrge Sprse Smmetrc Mtrces b the Precodtoed Cogte Grdet Method the Geerlzed Egele Problem 日 本 建 築 学 会 構 造 系 論 文 集 Vol.69 9-7 8.7 [4] 柏 木 光 博 ダブルシフト 逆 ベキ 乗 法 によるスパース 対 称 行 列 の 中 間 固 有 解 の 一 算 定 法 日 本 応 用 数 理 学 会 論 文 誌 Vol.9 No.33-38 9.9-3 -