旧バージョンの MATHEMATICAチュートリアル
最新の
Mathematica には,このチュートリアルに関連する新機能が追加されている.最新情報は
行列と線形代数を参照のこと.
線形代数の例題
このチュートリアルでは,線形代数を含むMathematica による計算の例を多数紹介する.
行列の並べ替え
疎行列技術によっては行列を並べ替え,要素をブロックにグループ化されるように試みる.その後各ブロックごとに密行列技術を使って計算が行われる.行列をブロックに並べ替える簡単な方法のひとつに,各行の要素の総和によってソートするというものがある.これを以下に示す.
| Out[6]= |  |
行列の置換を計算し,それを行列に適用する.これについては
「行列の置換」で述べた.
並べ替えられた行列は以下のようになる.図から並べ替えの効果が分かる.
| Out[9]= |  |
構造がなく,乱数の要素がたくさんある,はるかに大きい行列を考える.
| Out[17]= |  |
並べ替えられた行列のプロットから,要素がない行・列がたくさんあることが分かる.
| Out[18]= |  |
最大階数最小二乗解
行列方程式
を解くことは,行列計算の需要な問題のひとつである.その解法については「線形系の解法」で述べた.
が
×
の行列なら,
個の未知数を持つ
個の線型方程式一式である.
の場合は,方程式の数が未知数の数より多い過剰決定系である.最小二乗解を求める一般的な方法は「最小二乗法」で述べた.以下の例では,最大階数行列の過剰決定系の適切な解の別の求め方をいくつか紹介する.
このセクションの例題の解法は,入力として疎行列を与えたときにのみ正しく動作することに注意されたい.
最小二乗コレスキー
この方法は最小二乗解を求めるのに「コレスキー分解」を使う.行列
が最大階数を持つなら,解は以下のようにして求められる.
これらのステップは以下のようにして
Mathematica 関数に埋め込む.
| Out[4]= |  |
| Out[5]= |  |
| Out[6]= |  |
この方法は早いが,
を計算するため正確ではない.
求められた解の1-,2-,

ノルムを計算する.
| Out[7]= |  |
最小二乗QR
が最大階数の場合,問題をより正確に解くにはQR分解を使う.これを以下に示す.
解を求める
Mathematica プログラムは以下のようになる.
| Out[4]= |  |
QR分解の解はコレスキー法よりコストが高くなるが,より正確である.その差を以下に示す.
2-ノルムを使うとコレスキー解がどのくらい近くなるかを調べる.
| Out[9]= |  |
2-ノルムを使うとQR解がどのくらい近くなるかを調べる.QR分解の方がよりよい近似であることが分かる.
| Out[10]= |  |
1-ノルムと無限ノルムの最小化
行列方程式
が過剰決定形(
)の場合に,近似解を求める方法の多くでは,2-ノルムの最小化を行う.これは計算的に扱いやすくなる利点があるからである.ひとつの理由は,関数
が微分可能で,微分が線形操作であるということである.そのため,最小の解を求める線形系が形成できる.もうひとつの理由は,直交変換において2-ノルムが維持されるということである.これは,解くのがより簡単である可能性がある同様の問題が幅広く形成できるということを意味している.
しかし,1-ノルムや
-ノルム等,他のノルムを最小化することで見つかる他の解も存在する.これらは,解が個々の問題に関する重要な性質を維持する可能性があるので,特定の状況において,より望ましいことがある.次の例では,これらのノルムを最小化する近似解を求める方法を示す.両方とも制約付き線形問題の最小値を求める方法を使う.これは特に線形計画法として知られる.
Mathematica では,線形計画法の解法には関数LinearProgrammingが使える.これを使うと,整数,有理数,機械精度の近時実数,任意精度の近似実数等,Mathematica がサポートしている種々の型の数の線形計画問題が解ける.さらに,内点法を使って,特に大規模な系を最小化するのに適した方法も利用できる.
このセクションで与えられている解法は,入力が密行列の場合に適している.これを疎行列用に変更するのは簡単である.これは内点線形計画法を十分に利用するために必要である.
このセクションで述べる方法は,解の他の制約を加えて拡張することもできる.
1-ノルム最小化
1-を最小化する際は,次の式を最小化する
の値を求める.
これは新しい変数
を形成し,最小値を求めることで行える.
| Out[4]= |  |
求められた解の1-,2-,

ノルムを計算する.
| Out[5]= |  |
無限ノルムの最小化
-ノルムの最小化は1-ノルムの最小化と似ている.以下を最小化する
の値を求める.
これは新しい変数
を形成し,最小値を求めることで行える.
| Out[4]= |  |
求められた解の1-,2-,

-ノルムを計算する.
| Out[5]= |  |
有限差分解法
偏微分方程式を解くひとつの方法に,特別な離散化を行い,有限差分で解くというものがある.次の例で,単位正方形における基本的なポアソン(Poisson)方程式の差分解法を考える.

である80×80のグリッドの問題を設定する.変数行列が簡約された形で表示される.
Out[3]//Short= |
| |  |
ここで,変数行列の有限差分方程式の用途を作る.境界条件は境界付近のゼロ係数を持つ変数のひとつのオーバーハングを受け入れることで統合される.
Out[6]//Short= |
| |  |
ポアソン方程式の近似解は,これらの方程式を解くと求められる.方程式を解くひとつの方法は,
Solveを使うことである.
Solveは方程式を線形かつ疎なものであると認識し,疎行列ソルバを使用する.別のポアソン問題が解ける可能性のあるもうひとつの方法は,方程式から行列を作ることである.グリッドがベクトルに平坦化されると,行列は線形演算子を表す.
| Out[7]= |  |
| Out[8]= |  |
| Out[10]= |  |
同じ形状寸法で,右辺の関数は異なる場合に複数のポアソン方程式を解くには,行列の分解を計算し,それを繰返し使用する.
Mathematica で行列の分解を繰返し行う方法は
「行列の分解の保存」で述べた.
| Out[11]= |  |
グリッド上で関数
Cos[6x2+3y2]を表すベクトルを与える.
| Out[14]= |  |
| Out[17]= |  |
| Out[18]= |  |
| Out[19]= |  |
メッシュ区分
ここでは,固有値の実践的な使用方法を例示する.目的はメッシュ/グラフの頂点を2つのサブ領域に分割し,各部分の頂点の数をおおよそ同じにし,その境界が辺を切断するのをなるべく少なくするというものである(2つの頂点が別々のサブ領域にあると,辺は切断される).ラプラシアン(Laplacian)を作成し,そのフィードラー(Fiedler)ベクトル(2つ目に小さい固有値に相当する固有ベクトル)を使って頂点を並べるというものがある.これらはグラフ区分において,より効率的な方法であり,固有ベクトルの計算は必要ないが,フィードラーを使う方法は重要である.
データ
メッシュのプロット
メッシュは,各三角形をプロットするとプロットできる.
| Out[8]= |  |
ラプラシアン
グラフ
のラプラシアンは,
×
の疎行列である.
はグラフ中の頂点の数である.
の要素のほとんどはゼロであるが,対角要素は正の整数で,これは頂点の次数に対応している.ここで次数とは,ひとつの頂点につながっている他の頂点の数である.また,頂点
,
が辺でつながっている場合,要素
は-1である.
ここでのメッシュでは,非対角要素は疎配列
で与えられる.
| Out[10]= |  |
行の中の対角要素は非対角要素の和の正負を逆にしたものである.したがって,行列の対角要素は次のようになる.
| Out[11]= |  |
| Out[12]= |  |
| Out[13]= |  |
フィードラー(Fiedler)ベクトル
各行の要素の総和はゼロであるので,ラプラシアン行列は固有値がゼロの準正定値行列である.2つ目に小さい固有値に対応する固有ベクトルはフィードラーベクトルと呼ばれる.
固有値と固有ベクトルの最小のものから2つを求める.
| Out[15]= |  |
ノードの区分
フィードラー固有ベクトルを使ってベクトルの順序を求め,2つの頂点のグループ

と

を形成する.

の頂点の色を
Hue[0.](赤),

の頂点の色を
Hue[0.75](青)に設定する.
区分をプロットする.メッシュの特に密度の高い部分を避けるように区分され,各サブ領域はうまくつながっている.
| Out[23]= |  |
区分により辺が切断される回数は以下のようにして求められる.
| Out[26]= |  |
NDSolveを使った行列関数
Mathematica の数値ソルバの重要な機能のひとつに,入力がベクトルや行列でも使えるということが上げられる.一般にこれは連立方程式の解法に使われるが,行列特有の問題の解法にも使える.ここでは,数値微分方程式ソルバNDSolveを使って行列関数のパラメータかされた解を求める例を示す.
Out[2]//MatrixForm= |
| |  |
微分方程式の解は指数関数である.
| Out[3]= |  |
行列方程式を形成し
NDSolveで解くと,行列の指数
Exp[A t]が求められる.ここでの解は

のときにだけ有効である.
| Out[5]= |  |
解は
t についての行列関数
Exp[A t]を返す関数である.
t が0.0なら,単位行列が返される.
Out[6]//MatrixForm= |
| |  |
Out[8]//MatrixForm= |
| |  |
Mathematica には関数の指数を計算する関数
MatrixExpがある.これは求められた解と比較できる.
| Out[10]= |  |
| Out[11]= |  |
| Out[13]= |  |
Out[14]//MatrixForm= |
| |  |
| Out[16]= |  |
Out[17]//MatrixForm= |
| |  |

において恒等式

が成り立つことを示す.
Out[18]//MatrixForm= |
| |  |