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

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;6a=[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; endend, 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) • 線形方程式(系)

\ (バックスラッシュ) (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) Matlab任せ x=A\y; 10.8s 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)

最小二乗の各種解法 (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：逐次最小二乗に置き換え • 戦略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秒

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