Compose-3500：OMLでの線形代数演算

Tutorial Level: Beginner

ComposeのOML言語は、線形系の問題の解法、行列の解析、行列分解などを実行するための線形代数ツールボックスを提供します。

このチュートリアルでは、以下の操作方法について説明します：
  • 線形系を解く。
  • 基本的な行列解析を実行する。
  • 行列分解を実行する。

線形系の方程式の解法

  1. 線形方程式Ax = bの系をxについて解くには、\を直接使用するか、mldivide関数を使用します。どちらでも同じ結果が得られます。内部で選択されるアルゴリズムは行列Aに依存します。

    Aが正方行列の場合、OMLから一意な解が返されます。

    % solve linear equation A*x = b when A is square
% unique result is returned 
A = [9 13 5 2; 1 11 7 6; 3 7 4 1; 6 0 7 10];
b = [-1 6 3 -4]'; 
x1 = mldivide(A,b)
x2 = A\b
disp('check if x1 equals x2')
isequal(x1,x2)

    Composeでの結果を次に示します。

    Aに記述されている方程式の数が未知数の数より少ない場合、行が列よりも少なく、系は劣決定で、解が一意ではなくなります。計算結果は、次のように最小ノルムの最小二乗解になります。

    % solve under determined equation A*x = b
% there are fewer equations than unknowns in this case
% computed result is the least squares solution with the minimum norm
A = [1 1 1; 1 1 2];
b = [1 3]';
x1 = mldivide(A,b)
x2 = A\b
disp('check if x1 equals x2')
isequal(x1,x2)

    Composeでの結果を次に示します。

    Aに記述されている方程式の数が未知数の数より多い場合、行が列よりも多く、系は優決定で、解は存在しません。計算結果は、一意な最小二乗解になります。

    % solve over determined equation A*x = b
% A is rectangular and has more rows than columns
% the unique least squares solution is computed in this case
A = [2 1; -3 1; -1 1]; 
b = [-1 -2 1]'; 
x1 = mldivide(A,b)
x2 = A\b
disp('check if x1 equals x2')
isequal(x1,x2)

    Composeでの結果を次に示します。

    Aが非正則行列の場合は、次のように警告が表示されます。

    % solve linear equation A*x = b when A is singular
% it will raise a warning in this case
A = [16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1];
b = [1 1 1 1]';
x1 = mldivide(A,b)
x2 = A\b
disp('check if x1 equals x2')
isequal(x1,x2)

    Composeでの結果を次に示します。

  2. Ax = bを解く際、行列Aに関する演繹的知識（下三角行列、正定値行列、一般矩形行列など）がある場合は、linsolve関数を選択することで効率に優れた演算を実現できます。基本的に、linsolveは\と同様に機能しますが、行列Aの特性を指定する入力引数があります。
  3. ここではAx = bを解いていますが、xA = bを解く場合もあります。この場合は、前の項と同じロジックで/またはmrdivide関数を使用できます。

基本的な行列解析

本項では、行列解析の2つの事例を紹介します。1つは逆行列、もう1つは固有値と固有ベクトルです。

  1. 逆行列を計算します。
    行列が正方で正則の場合は、inv関数を使用して逆行列を計算します。
    % calculate non-singular matrix inversion
format(3,2);
M = [1 2 -3;4 5 6;-7 8 0]
Inv_M = inv(M)
M2 = inv(inv(M))
find((M - M2) > 10*eps)

    Composeでの結果を次に示します。

    上の図で、inv(inv(M))が、丸め誤差の範囲内でMと等しいことがわかります。なお、Mが正方行列ではない場合または非正則行列である場合は、invからエラーまたは警告が返されます。それ以外の場合は、そのムーア・ペンローズ逆行列または一般逆行列を計算できます。A2がA1の一般逆行列である場合は、次の条件を満たす必要があります。

    A2*A1*A2 = A2
A1*A2*A1 = A1
conj(A1*A2) = A1*A2 
conj(A2*A1) = A2*A1

    OMLの4つの規則すべてに従っていることを確認します。

    % calculate singular matrix inversion using pinv
format(3,2)
A1 = [1 2 3;4 5 6;7 8 9];
disp('A2 vs. A2*A1*A2')
A2 = pinv(A1)
A2*A1*A2
disp('A1 vs. A1*A2*A1')
A1
A1*A2*A1
isequal(conj(A1*A2),A1*A2)
isequal(conj(A2*A1),A2*A1)

    Composeでの結果を次に示します。

  2. 行列の固有値と固有ベクトルを計算します。
    正方行列の固有値は、その行列による線形変換を経ても方向が変化しないベクトルです。これは次の式で表すことができます。

    は、固有ベクトルに関連付けられた固有値であるスカラーです。OMLでは、eig関数を使用して固有値を計算します。

    % calculate eigenvalue and eigenvector of a square matrix
A = [2 0 0;0 3 4;0 4 9];
[V,D] = eig(A)
% verify result satisfy A*V = V*D
find((A*V - V*D) > eps)

    Composeでの結果を次に示します。

行列分解

行列分解は、行列を行列の積に因数分解する演算です。さまざまな行列分解がありますが、それぞれが特定のクラスでの問題を対象としています。本項ではLU分解とQR分解を紹介します。

  1. LU分解

    LU分解は、行列を下三角行列と上三角行列の積として因数分解します。OMLでは、次のようにlu関数を使用します。

    % LU decomposition
A = [4 3 1;0 6 3;1 2 4]
[L,U] = lu(A)
% verify if 
L*U = A

    Composeでの結果を次に示します。

    この関数から返されたLとUがそれぞれ下三角行列と上三角行列で、L*Uが元の行列Aと等しいことがわかります。

  2. QR分解

    QR分解は、行列を直交行列Qと上三角行列Rの積として因数分解します。OMLでは、qrを使用します。

    % QR decomposition
format short
A = [12 -51 4; 6 167 -68; -4 24 -41]
[Q,R] = qr(A)
% verify if Q is an orthogonal matrix
assert((Q'*Q - eye(size(A))) < 10*eps)
assert((Q*Q' - eye(size(A))) < 10*eps)
% verify if Q*R = A
assert((A - Q*R) < 100*eps)

    Composeでの結果を次に示します。

    Qが直交行列、Rが上三角行列で、Q*Rが丸め誤差の範囲内でAと等しいことがわかります。

  3. スクリプト全体については、install_dir/tutorials/をご参照ください。