octave 擬合線性方程式範例改用 SVD

% 假設實驗數據 (x̂,ŷ) 計算最小誤差線性方程式 ŷ = c1 + c2*x̂ + c3*x̂² + c4*x̂³ 的解

samples = 64; % 測量樣本數

c1 = 10;

c2 = 0.8;

c3 = 0.5;

c4 = 0.3;

f = @(t) c1 + c2*t + c3*(t^2) + c4*(t^3);

y = [];

x1 = [];

x2 = [];

x3 = [];

for (t = 1 : samples) % 取樣, Δt

s = f(t);

err = rand - 0.5; % 誤差訊號

noise = s * err; % max +- 50% signal

y = [y, s + noise]; % 讀取測量訊號組成向量

x1 = [x1, t];

x2 = [x2, t^2];

x3 = [x3, t^3];

end

elements = 4

P = [ones(1, samples); x1; x2; x3]';% 測量矩陣

% pseudo inverse to solve Ax=b => x = A⁺b

[U,S,V] = svd(P); % A = U * Σ * V'

sigma = S*0; % it's not possible to inv(S), so copy full matrix first

for (t = 1 : elements)

sigma(t,t) = 1/S(t, t); % S⁻ is just reciprocal, fill at diagonal position

end

R = V * sigma' * U'; % R = A⁺ = V * Σ⁻ * U'

M = R * y';

g = []

for (t = 1 : samples)

z = 0;

for (el = 1 : elements)

z = z + M(el) * (t^(el-1));

end

g = [g, z];

end

plot(x1, y, 'g:x', x1, g(x1) ,'r-');

備註:

Jacobi SVD 演算法參考

https://netlib.org/lapack/lawnspdf/lawn15.pdf

在第 32 頁 Algorithm 4.1 One-sided Jacobi for the singular value

針對演算法, A 必需是 nxn 對稱的方陣, 我的理解是:
1. G = A 原矩陣, V = I 單位方陣(對角線都是 1, 其他地方為 0, 行列與 A 相同都是 nxn 方陣)
2. 計算每一對行向量 (Gî, Gĵ), 分量分別是 Gki, Gkj, 其中 k = 0 : n-1
a = Gî•Gî = Σ Gim² = |Gî|² => 行向量 i 的自內積 = 行向量 i 的長度平方
b = Gĵ•Gĵ = Σ Gjm² = |Gĵ|² => 行向量 j 的自內積 = 行向量 j 的長度平方
c = Gî•Gĵ = Σ(Gim*Gjm) => 兩行向量 i,j 的內積
ζ = (b - a) / (2 * c) => 有可能是負值
t = sign(ζ) / (|ζ| + √(1+ζ²)) => 類似 tanθ, 直角三角形, 邊長分別是: t,1,√(1+t²)
cs = 1 / √(1 + t²) => 類似 cosθ, 餘旋邊長 = 1, 斜邊長=√(1+t²)
sn = cs * t => 類似 sinθ, 正旋邊長 = t, tanθ = sinθ/cosθ = t
T = [cs, -sn;sn cs] => 就是一個旋轉矩陣
3. 透過旋轉矩陣 T 與向量相乘, 同時更新行向量組 (Gî, Gĵ), (Vî, Vĵ), 也就是把座標軸加以旋轉:
(G'î, G'ĵ) = T (Gî, Gĵ)
(V'î, V'ĵ) = T (Vî, Vĵ)
|Gkî,Gkĵ|ᵗ 疊代更新似乎是 2x2 旋轉矩陣 與 前一次 |Gkî,Gkĵ|ᵗ 行向量內積的結果:
|G[k][i]| = (cs , -sn) • |G[k][i]|
|G[k][j]| = (sn , cs) • |G[k][j]|
|Vkî,Vkĵ|ᵗ 疊代更新似乎是 2x2 旋轉矩陣 與 前一次 |Vkî,Vkĵ|ᵗ 行向量內積的結果:
|V[k][i]| = (cs , -sn) • |V[k][i]|
|V[k][j]| = (sn , cs) • |V[k][j]|
==============================================
4. 更新直到 |c|/√(ab) <= tol 可接受的值
|c| = |Gi•Gj|/(|Gi|²|Gj|²) = |cosθ|
當所有夾角接近 90 度時, cosθ 逼近於 0, 就能結束, 也就是行向量將會彼此互相垂直 ?