宮崎大学 >> 工学部 >> 情報システム工学科 >> 伊達 >> 講義で使う Octave >> 線形代数 >>
Last modified: Mon Nov 8 11:23:15 JST 2010

Octave と線形代数

目的

講義「線形代数」では，日常生活では使わない，難しそうな用語を使う．しかし、用語に騙されてはいけない。それぞれの用語が意味することは、実は、単純だ。これが本当だということを実感してほしい。

コンピュータを操作するだけで，以下の用語の概念、概念どうしの関係に対する、直観的な理解をつかもう。

行列式
逆行列
階数（ランク）
固有値・固有ベクトル

どのくらいこの方法が有効か、あとで感想を聞かせてください。
もちろん教科書を何回も読んで、それぞれの概念を数学的に理解する必要はあります。

ここでは行列とベクトルの演算が楽にできる octave とよばれるソフトウェアを使ってみる

まず linux を立ち上げ，octave を起動します．

% octave

GNU Octave, version 3.0.5
Copyright (C) 2008 John W. Eaton and others.
This is free software; see the source code for copying conditions.
There is ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  For details, type `warranty'.
.
.
For information about changes from previous versions, type `news'.

octave:1>

電卓と同じように使えることを確認しよう。例: 1+2 （=は入力しない）
```
octave:2> 1+2
ans =  3
```

2x2 の行列 A を定義する．

octave:3>  A=[1.0, -0.3;  -0.7, 0.6]

A =

   1.00000  -0.30000
  -0.70000   0.60000

2次元ベクトル x を定義する．

octave:4>  x = [-3.0, 0.5]

x =

  -3.00000   0.50000

行列とベクトルの積、Ax を求めてみよう。 Ax を手で計算すると、( 1*(-3) + (-0.3)*0.5, -0.7*(-3)+0.6*0.5 ) = (-3.15, 2.4)
```
octave:5>  A*x
```
```
error: operator *: nonconformant arguments (op1 is 2x2, op2 is 1x2)
error: evaluating binary operator `*' near line 8, column 2
```
こうすると怒られる．x は横ベクトルであるので、積 Ax が定義できない。こういう理由で怒られた．
横ベクトル x を転置して、縦ベクトルにしておく必要がある。
x' で、x の転置を作ることができる。（アポストロフィが転置を意味する．）まず x を転置して，縦ベクトルにしておく．
```
octave:6>  x = x'
```
```
x =

  -3.00000
   0.50000
```
こうすることで、確かに x が縦ベクトルになった．
Ax を求める．
```
octave:7>  A*x
```
```
ans =

  -3.1500
   2.4000
```
手で計算した結果と一致している。

行列式，逆行列、ランク（階数）、固有値
以下ではすべて、正方行列 A について話しを進める。
行列 A の行列式は det(A)，逆行列は inv(A)，ランク（階数）は rank(A)，固有値は eig(A) で計算できる！
（Aを A=[1.0, -0.3; -0.7, 0.6] などと具体的に与えて、キーボードを叩くだけ．）
単に A と入力すれば，現在、変数Aに代入されている値（行列）が表示される．
```
octave:8>  A
```
```
A =

   1.00000  -0.30000
  -0.70000   0.60000
```
det(A) と入力すると、行列Aの行列式の値が返ってくる。
```
octave:9> det(A)
ans =  0.39000
```
同様に inv(A) と入力すると、行列Aの逆行列が返ってくる。
```
octave:10> inv(A)
ans =

   1.53846   0.76923
   1.79487   2.56410
```
A の逆行列を B としておこう。 AB = BA = I である。（これは定義。I は単位行列）
```
octave:11> B = inv(A)
B =

   1.53846   0.76923
   1.79487   2.56410
```
ここで、「行列は写像だ」、「行列式は、面積拡大率・体積拡大率」であることを思いだす。
BAx という変換（BはAの逆行列）は、x をAで変換して、その後、もとに戻している（BAx = Ix）。
もとに戻したのだから、面積拡大率は1である。つまり |B||A|=1になる。それを確認しよう。
```
octave:12> det(B)
ans =  2.5641
```
```
octave:13> det(A)*det(B)
ans =  1
```
Aのランク（階数）を、もとめてみよう。
```
octave:14> rank(A)
ans =  2
```
この行列 Aは、筋のいい行列。筋の悪いのも試しておこう。
A = [ 1, 2; 2, 4] なら、rannk(A) は、どうなる？
A = [ 0, 0; 0, 0] なら、rannk(A) は、どうなる？
Aの固有値を、もとめてみよう。eig という関数を使う。
```
octave:15> eig(A)
ans =

   1.30000
   0.30000
```
これは、行列Aの固有値は2つあって、1.3 と 0.3 ということを意味している。
数学的に勉強した人にはあたりまえだが、実は、固有値は、固有ベクトルとペアになっている。つまり、行列Aには、固有値1.3 に対応する固有ベクトルが1つ、固有値0.3 に対応する固有ベクトルが1つ存在する。これをみてみよう。
eig という関数を用い固有ベクトルを求める．
```
octave:16>  [V, LAMBDA] = eig (A)
```
```
V =

   0.70711   0.39392
  -0.70711   0.91915

LAMBDA =

Diagonal Matrix

   1.30000         0
         0   0.30000
```
行列 V は固有ベクトル（縦ベクトル）を2つ並べたものになっている。教科書の問題を解くと、固有ベクトルの要素は簡単な整数で表現できるようになっているが、コンピュータはそんなことは知らない。実は、ここで表示された固有ベクトルは、わざわざ大きさが1に正規化されている．（固有ベクトルの長さは、重要ではないこということ！）
さて、行列 V の第一列（固有ベクトル1）だけをとりだし，それを V1 とする．
```
octave:17> V1=V(:,1)
```
```
V1 =

  -0.70711
  -0.70711
```
ベクトル V1 の大きさ（ノルム、長さを一般化したもの）を計算する．
```
octave:18> norm(V1)
```
```
ans = 1
```

3x3の行列 A を定義する．

octave:19>  A=[5,-1,-1; 1,2,0; 3,-1,1]

eig という関数を用い固有ベクトルを求める．

octave:20>  [V, LAMBDA] = eig (A)

V =

   -5.7735e-01    5.7735e-01   -1.9025e-15
   -5.7735e-01    5.7735e-01   -7.0711e-01
   -5.7735e-01    5.7735e-01    7.0711e-01

LAMBDA =

  3.00000  0.00000  0.00000
  0.00000  3.00000  0.00000
  0.00000  0.00000  2.00000

【補足】この -5.7735e-01 という数は何か．

octave: 21> 1.2345*0.00000000000000000001
ans =  1.2345e-20

としてみれば何を意味しているかわかる．（10 の -20乗ということ）

octave には、いろいろなコマンド、関数がある。 google で「octave 使い方」とすれば、だいたいのことは解決する。

もちろんオンラインマニュアルもある。例えば、 eig という関数の使い方について知りたければ

octave:22> help eig

eig is the dynamically-linked function from the file
/usr/libexec/octave/2.1.73/oct/i686-pc-linux-gnu/eig.oct

 -- Loadable Function: LAMBDA = eig (A)
 -- Loadable Function: [V, LAMBDA] = eig (A)
     The eigenvalues (and eigenvectors) of a matrix are computed in a
     several step process which begins with a Hessenberg decomposition,
     followed by a Schur decomposition, from which the eigenvalues are
     apparent.  The eigenvectors, when desired, are computed by further
     manipulations of the Schur decomposition.

英語を読もう。スペースキーで続きが読める。q でおしまい。

行列 A の 100乗を計算する．

octave: >  A=[3,2; 4,1]

octave: >  A^100

ans =

   5.2591e+69   2.6295e+69
   5.2591e+69   2.6295e+69

100 x 100 の行列を作成してみる．各要素にはでたらめな値を入力しておく。 randn は正規分布乱数を成分にする行列を生成．以下では、最後に ; を付けないと生成した行列（大量の数字）を表示してしまう．止める場合には q キーを押す．
```
octave: >  M = randn(100,100);
```
```
octave: >  [P, LAMBDA] = eig (M);
```
この程度の大きさの行列の固有値なら簡単に計算できる！
そのほか1： 3x3 の単位行列 I の作りかた。
```
octave: >  I = eye(3,3)
```
そのほか2：対角行列の作りかた。
```
octave: >  A=diag([1, 2, 3])
```
ほかにも、zeros(3,3)、ones(3,4) などで、単純な構造をもつ行列が生成できる。

とりあえず、ここまで。あとは余裕があれば読んでみる。　

計算機を使うと，対角化などとは関係なく，簡単に結果が得られるが...

octave: >  [P, LAMBDA] = eig (A)

P =

   0.70711  -0.44721
   0.70711   0.89443

LAMBDA =

   5   0
   0  -1

octave: > D = inv(P)*A*P

D =

    5.0000e+00    3.4022e-16
   -1.1471e-16   -1.0000e-00

非対角成分が 0 にはなってない（限りなく0に近い，小さな値ではあるが）．このあたりの事情は数値計算の講義で学ぶ（と思う）．

octave: > P*(D^100)*inv(P)

ans =

   5.2591e+69   2.6295e+69
   5.2591e+69   2.6295e+69

octave: >  A^100

ans =

   5.2591e+69   2.6295e+69
   5.2591e+69   2.6295e+69

めでたし？

リンク

宮崎大学 >> 工学部 >> 情報システム工学科 >> 伊達 >> 講義で使う Octave >>

Octave と線形代数

目的

ここでは行列とベクトルの演算が楽にできる octave とよばれるソフトウェアを使ってみる

行列式，逆行列、ランク（階数）、固有値

リンク