Transcript Document

MATLAB測位プログラミングの
基礎とGT (3)
東京海洋大学産学官連携研究員
高須 知二
行列演算プログラミング
•
•
•
•
•
行列生成
行列操作
行列演算
\ (バックスラッシュ)
組込関数
行列生成 (1)
• 空 : a=[];
• スカラー (0次元配列) :
a=1; b=3+3i; c=pi; d=NaN; c=Inf;
• ベクトル (1次元配列) : a=[1 2]; b=[3;4;5]
• 行列 (2次元配列) : a=[1 2 3 4; 5 6 7 8];
• (3階以上) テンソル(3次元以上配列)
a=cat(3,[1 2 3;4 5 6],[7 8 9;10 11 12]);
行列生成 (2)
•
•
•
•
零行列 : a=zeros(10); a=zeros(4,5);
単位行列 : a=eye(10);
等間隔(行)ベクトル : a=1:10; b=0:-0.5:-10;
単一値行列:
a=ones(10); a=4*ones(10);
a=repmat(NaN,10,10); a(1:10,1:10)=Inf;
• 乱数 : a=rand(10); a=randn(10,1);
行列生成 (3)
• 要素への直接代入
a(1,1)=1; a(1,2)=2; a(1,3:7)=3;
• 行列サイズ自動拡張 (ゼロfill)
clear; a(10,8:10)=1;
• 行列の連結
a=[1 2 3]; b=[4 5]; c=[a b];
d=[1:4 10:12 8:-1:6];
e=[[1;2;3];[[4;5;6],[7;8;9]]];
行列操作 (1)
• 要素参照・代入
a=[1 2 3 4;5 6 7 8;9 10 11 12];
b=a(1,3); c=a(:,1); d=a(2,:); e=a(1:3,1:2);
f=a([1 2],[1 3]); g=a(end,2); h=(1,end-2);
k=a([1 1 1 1 1],[1 2 1 2]);
a(1,1)=-1; a(1,2:3)=0; a(end,:)=3:6;
• 行列サイズ、次元
[n,m]=size(a); n=length(b); n=ndims(c);
行列操作 (2)
• 転置
a=[1 2 3;4 5 6;7 8 9]; b=a'; c=(1:10)'; (or .')
• サイズ変更
b=reshape(a,9,1);
c=a(:); d=zeros(1,9); d(:)=a;
• 交換
b=a(:,[2 1 3]); c=a([3 2 1],[3 1 2]);
d=flipud(a); e=fliplr(a);
行列操作 (3)
• 要素追加
a=[1 2 3;4 5 6;7 8 9];
b=[a;[10 11 12]]; b=[10;11;12;a];
a(end+1,:)=[10 11 12];
• 要素削除
a(:,1)=[]; a(end,:)=[];
行列操作 (4)
• FIND() : 非ゼロ要素のインデックス
a=[0 0 0 1 0 1]; i=find(a); i->4;6
a=[2 0;0 1;0 0]; [i,j]=find(a); i->1;2, j=1;2
• ベクトル演算関数、演算子との組み合わせ
→forループ代替 : matlab表現, 実行効率
• 行列要素参照中ではfind()を省略可能
a(find(a<0)) <-> a(a<0)
a(find(a(:,1)==3),end) <-> a(a(:,1)==3,end)
行列操作 (5)
• 例1 : 0.5未満の要素の個数を数える
(1) n=0; for i=1:length(a), if a(i)<0.5, n=n+1; end, end, n
(2) n=length(find(a<0.5))
(3) n=sum(a<0.5)
• 例2 : NaN以外の値の平均を求める
(1) s=0;n=0;
for i=1:length(a)
if ~isnan(a(i)), s=s+a(i); n=n+1; end
end, m=s/n
(2) m=mean(a(~isnan(a)))
行列操作 (6)
• 例3 : for文とfindの組み合わせ
a=rand(100,4);
for i=find(a(:,1)<0.5)'
disp(sprintf('%f',a(i,:)))
end
行列演算 (1)
•
•
•
•
•
•
•
加減算：行列+スカラー : A+2, 3+A
加減算：行列+行列 : A+B
乗算：行列×スカラー : A*4, 5*A
乗算：行列×行列 : A*B
要素乗算：行列×行列 : A.*B
べき乗 : A^3 =A*A*A (A:正方行列)
要素べき乗：A.^3
行列演算 (2)
• 比較 (==, ~=, <, >, <= >=)
行列 : 行列
行列 : スカラー
結果は(0,1)要素行列で得られる→FIND()
• 零行列判定 : isempty(a) (×a==[])
• 行列内容一括比較
c=isequal(a,b) (サイズ+内容)
d=all(all(a>b)); e=any(any(a<b));
\ (バックスラッシュ) (1)
• 線形方程式(系)
ym1  Amn xn1
y1  a11x1  a12x2  a13x3  ...  a1n xn
y2  a21x1  a22x2  a23x3  ...  a2n xn
y3  a31x1  a32x2  a33x3  ...  a3n xn
...
ym  am1 x1  am2 x2  am3 x3  ...  am nxn
\ (バックスラッシュ) (2)
• matlabでの解法 : x=A\y
(1) n=m : 線形方程式→ ガウス消去法 等
(2) n<m : 最小二乗解→QR分解
(3) n>m : 一般解の一つ (basic解?)
• Aが三角、対称、スパースか否か等を判定
し最適(精度、効率)な解法を内部選択
• x=y/A (スラッシュ)→y=x*Aの解
A/B = (B'\A')'
\ (バックスラッシュ) (3)
• 例1 : 観測値の2次多項式フィッティング
(1) A=[x.^2 x ones(size(x))]; a=A\y;
(2) a=polyfit(x,y,2);
• 例2 : 単独測位アルゴリズム
最小二乗各種解法 (1)
x=A\y;
x=(A*A')\(A'*y);
x=inv(A*A')*(A'*y);
R=chol(A‘*A); x=R\(R'\(A'*y));
[Q,R]=qr(A,0); x=R\(Q'*y);
10.8s Matlab任せ
5.3s 正規方程式
6.0s
正規方程式
5.2s コレスキ分解
11.8s
QR分解
[U,D,V]=svd(A,0); x=V*(D\(U'*y)); 46.4s 特異値分解
x=pinv(A)*y;
51.9s 疑似逆行列
(n=1000, m=5000, Pentium 4 3.2GHz, Matlab 6.5.1)
最小二乗の各種解法 (2)
• Aがランク落ちしていた場合の動作の違い
(1) x=A\y;
→ ランク落ち警告+ basic解
(2) x=inv(A'*A)*(A'*y)
→ エラーまたは不正解
(3) x=pinv(A)*y;
→警告無し、ノルム最小解
• 基本的に(2)は使うべきでない。
大規模最小二乗問題 (1)
•
•
•
•
パラメータ数>数1000
観測データ数>数100000
計画行列が実メモリ上に載り切らない
実用的なパラメータ推定問題ですぐ現れる
→
計算は結構やっかい
大規模最小二乗問題 (2)
• 戦略1：逐次最小二乗に置き換え
1
x  ( A A) A y  x 
T
T
( A A ) ( A
1
T
i
i
i
T
yi )
• 戦略2 : カルマンフィルタに置き換え
• 戦略3 : スパース計画行列 + matlab + \
• 戦略4 : 最小二乗演算ライブラリを使う
行列用組込関数 (1)
• 行列用演算 :
diag(), norm(), rank(), det(), trace(), inv(),
dot(), cross(), meshgrid()... (see help)
• 行列入力→行列出力関数 :
sin(), cos(), tan(), sqrt(), exp(), log(), floor(),
mod(), ...
(ほとんどの初等関数、特殊関数)
• スカラ計算式→行列計算式
行列用組込関数 (2)
• 例1 :
• 例2 : 2次元メッシュ/3Dグラフ生成
[x,y]=meshgrid(0:0.1:20,0:0.1:20);
surf(x,y,sin(x)+cos(y))
行列演算高速化 (1)
• 行列サイズの事前割当
→サイズ自動拡張をなるべく使わない
• 行列のまま計算、find()の利用
→forループをなるべく使わない
• find()の高速化 :
インデックス操作の工夫
sort(), sortrows(); unique(), intersect()
行列演算高速化 (2)
• 例1 : 100万回ループ
x=0:0.1:100000; y=zeros(1000000,1);
y=sin(x);→0.4秒
for i=1:length(x), y(i)=sin(x(i)); end;→3秒
• 例2 : 100万回ループ
x=zeros(1000000,1);
for i=1:length(x), x(i)=i; end→1.4秒
x=[]; for i=1:1000000, x=[x;i]; end>1000秒