MATLAB測位プログラミングの 基礎とGT (3)

Presentation on theme: "MATLAB測位プログラミングの 基礎とGT (3)"— Presentation transcript:

1 MATLAB測位プログラミングの 基礎とGT (3)
東京海洋大学産学官連携研究員 　高須　知二

2 行列演算プログラミング 行列生成 行列操作 行列演算 \ (バックスラッシュ) 組込関数

3 行列生成 (1) 空 : a=[]; スカラー (0次元配列) : a=1; b=3+3i; c=pi; d=NaN; c=Inf;
(3階以上) テンソル(3次元以上配列) a=cat(3,[1 2 3;4 5 6],[7 8 9; ]);

4 行列生成 (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);

5 行列生成 (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]]];

6 行列操作 (1) 要素参照・代入 a=[ ; ; ]; 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([ ],[ ]); a(1,1)=-1; a(1,2:3)=0; a(end,:)=3:6; 行列サイズ、次元 [n,m]=size(a); n=length(b); n=ndims(c);

7 行列操作 (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);

8 行列操作 (3) 要素追加 a=[1 2 3;4 5 6;7 8 9]; b=[a;[ ]]; b=[10;11;12;a]; a(end+1,:)=[ ]; 要素削除 a(:,1)=[]; a(end,:)=[];

9 行列操作 (4) FIND() : 非ゼロ要素のインデックス a=[ ]; 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)

10 行列操作 (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)))

11 行列操作 (6) 例3 : for文とfindの組み合わせ a=rand(100,4); for i=find(a(:,1)<0.5)' disp(sprintf('%f',a(i,:))) end

12 行列演算 (1) 加減算：行列+スカラー : A+2, 3+A 加減算：行列+行列 : A+B 乗算：行列×スカラー : A*4, 5*A
べき乗 : A^3 =A*A*A (A:正方行列) 要素べき乗：A.^3

13 行列演算 (2) 比較 (==, ~=, <, >, <= >=) 行列 : 行列 行列 : スカラー 結果は(0,1)要素行列で得られる→FIND() 零行列判定 : isempty(a) (×a==[]) 行列内容一括比較 c=isequal(a,b) (サイズ+内容) d=all(all(a>b)); e=any(any(a<b));

14 \ (バックスラッシュ) (1) 線形方程式(系)

15 \ (バックスラッシュ) (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')'

16 \ (バックスラッシュ) (3) 例1 : 観測値の2次多項式フィッティング (1) A=[x.^2 x ones(size(x))]; a=A\y; (2) a=polyfit(x,y,2); 例2 : 単独測位アルゴリズム

17 最小二乗各種解法 (1) 10.8s Matlab任せ 5.3s 正規方程式 6.0s 正規方程式 5.2s コレスキ分解 11.8s
x=A\y; 10.8s Matlab任せ x=(A*A')\(A'*y); 5.3s 正規方程式 x=inv(A*A')*(A'*y); 6.0s 正規方程式 R=chol(A‘*A); x=R\(R'\(A'*y)); 5.2s コレスキ分解 [Q,R]=qr(A,0); x=R\(Q'*y); 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)

18 最小二乗の各種解法 (2) Aがランク落ちしていた場合の動作の違い (1) x=A\y; 　　→ ランク落ち警告+ basic解 (2) x=inv(A'*A)*(A'*y) 　　→ エラーまたは不正解 (3) x=pinv(A)*y; 　→警告無し、ノルム最小解 基本的に(2)は使うべきでない。

19 大規模最小二乗問題 (1) パラメータ数>数1000 観測データ数>数100000 計画行列が実メモリ上に載り切らない
実用的なパラメータ推定問題ですぐ現れる → 計算は結構やっかい

20 大規模最小二乗問題 (2) 戦略1：逐次最小二乗に置き換え 戦略2 : カルマンフィルタに置き換え
戦略3 : スパース計画行列 + matlab + \ 戦略4 : 最小二乗演算ライブラリを使う

21 行列用組込関数 (1) 行列用演算 : diag(), norm(), rank(), det(), trace(), inv(), dot(), cross(), meshgrid()... (see help) 行列入力→行列出力関数 : sin(), cos(), tan(), sqrt(), exp(), log(), floor(), mod(), ... (ほとんどの初等関数、特殊関数) スカラ計算式→行列計算式

22 行列用組込関数 (2) 例1 : 例2 : 2次元メッシュ/3Dグラフ生成 [x,y]=meshgrid(0:0.1:20,0:0.1:20); surf(x,y,sin(x)+cos(y))

23 行列演算高速化 (1) 行列サイズの事前割当 →サイズ自動拡張をなるべく使わない
行列のまま計算、find()の利用 →forループをなるべく使わない find()の高速化 : インデックス操作の工夫 sort(), sortrows(); unique(), intersect()

24 行列演算高速化 (2) 例1 : 100万回ループ x=0:0.1:100000; y=zeros( ,1); y=sin(x);→0.4秒 for i=1:length(x), y(i)=sin(x(i)); end;→3秒 例2 : 100万回ループ x=zeros( ,1); for i=1:length(x), x(i)=i; end→1.4秒 x=[]; for i=1: , x=[x;i]; end>1000秒