延續上一篇利用 Givens 旋轉矩陣 Gⱼ(θ), 將 GⱼᵗG 不斷夾在 IA之間, 產生兩個新矩陣, 最後在右邊產生上三角矩陣 R, 左邊自然生成一個正交矩陣 Q 的過程
A = IA = I(GⱼᵗGⱼ)A = (IGⱼᵗ ...) * (... Gⱼ)A = Q * R
詳細參考 https://masontseng.blogspot.com/2023/04/c-givens-qr.html
Jacobi 旋轉矩陣 Jⱼᵢ(θ) 類似 Givens 旋轉矩陣 Gⱼ(θ), 也是在單位矩陣的對角線上,將 2x2 旋轉矩陣的元素擺放到 4
個角落: J(j,j)=cosθ, J(i,j)=-sinθ, J(j,i)=sinθ, J(i,i)=cosθ, 而且Jacobi 矩陣的轉置矩陣也是反矩陣 Jᵗ(θ)=J(θ)⁻, 也就是說 J(θ)Jᵗ(θ) = I. 利用 Jacobi 旋轉矩陣不停將一對 JⱼᵢJⱼᵢᵗ 夾在 AI 中間,產生兩個新矩陣, 詳細分解過程如下:
A = AI = A(J₁J₁ᵗ)I = AJ₁ * J₁ᵗI = (AJ₁) * (IJ₁)ᵗ
= (AJ₁)J₂J₂ᵗ(IJ₁)ᵗ = (AJ₁J₂) * (IJ₁J₂)ᵗ
= (AJ₁J₂) ... JⱼᵢJⱼᵢᵗ ... (IJ₁J₂)ᵗ
= (AJ₁J₂J₃ ... Jⱼᵢ) * (IJ₁J₂ ...Jⱼᵢ)ᵗ
= L * Vᵗ = (US)Vᵗ = USVᵗ
由此可知 A 和 I 只要分別不斷往右同時乘上一個旋轉矩陣 Jⱼᵢ(θ), 就可同時生出兩個新矩陣 L 和 Vᵗ, 這過程稱為單邊 Jacobi SVD (Singular Value Decomoposition) 分解法, 因Jacobi 矩陣 Jⱼᵢ(θ) 是正交矩陣, 透過不斷插入一對 JⱼᵢJⱼᵢᵗ 新矩陣, 右邊自然分解出一個正交矩陣 Vᵗ, 結束後, 左邊 L 矩陣再透過對角化,分解成兩個新矩陣 U 和 S(只要將 U 變成單位向量所組成的矩陣, S 就是行向量長度所組成的對角線矩陣), 為了要把左邊矩陣 L 化成正交矩陣, 旋轉角計算方式就要特別設計, 目的是將所有 L = [L̂] 內的行向量 L̂ₖ, 一一旋轉變成相互正交的行向量 L̂ᵢ, L̂ⱼ, 也就是讓 L̂ᵢ•L̂ⱼ = 0 , 之後算出每個行向量 L̂ₖ 的長度 σₖ =||L̂ₖ||₂ , 讓所有正交的行向量歸一化, 變成單位行向量 ûₖ = L̂ₖ/σₖ (從這不難看出 σ 稱為奇異值 singular value 的原因, 若 σₖ 等於 0, ûₖ 將變成無窮大的奇向量, 因此必須把奇異值 σ = 0 這一項剔除掉, 對應的 ûₖ/v̂ₖ 行向量也順勢移除, 或者將他排序過, 將奇異值 0 擺至矩陣最後面, 讓它不影響計算結果), 最後將正交單位行向量組成 U = [û] 矩陣, 奇異值 σₖ 擺至另一個矩陣的對角線上, 組成 S 矩陣: σₖ = ||L̂ₖ||₂, S = Λ = [σ].
任何矩陣都可以透過 Jacobi SVD分解成 3 個 U S V 矩陣, 演算法詳見以下在第 32 頁 Algorithm 4.1 One-sided Jacobi for the singular value 文章 https://netlib.org/lapack/lawnspdf/lawn15.pdf . 底下是 Jacobi SVD 用疊代法分解的程序,算出 L, V 矩陣:
初始化 L = A, V = I, Jⱼᵢ(θ)=[cosθ, -sinθ ; sinθ, cosθ]
L = L *Jⱼᵢ(θ)
V = V *Jⱼᵢ(θ)
// svdTest.c
#include <stdio.h>
#include <math.h>
#include <vector>
#define tol ((double)1e-12)
using namespace std;
using F_vector = vector<double>;
using F_matrix = vector<F_vector>;
void vectorShow(F_vector v) {// 列出向量元素
for (int i = 0; i < v.size(); i ++) { printf("%-10.4f \t", v[i]); }
printf("\n");
}
void matrixShow(F_matrix matrix) {// 列出矩陣元素
for (int j = 0; j < matrix.size(); j ++) {
printf("<%-2d:>\t", j);
vectorShow(matrix[j]);
}
printf("\n");
}
F_matrix transpose(F_matrix &matrix) { // 矩陣轉置
int row = matrix.size();
int col = matrix[0].size();
F_matrix temp(col, F_vector(row));
for (int i = 0; i < row; i++) {
for (int j = 0; j < col; j ++) {
temp[j][i] = matrix[i][j];
}
}
return temp;
}
F_matrix eye(int n) {// 單位矩陣
F_matrix temp(n, F_vector(n));
for (int i = 0; i < n; i ++) {
for (int j = 0; j < n; j++) {
temp[i][j] = (i == j);
}
}
return temp;
}
F_matrix product(F_matrix q, F_matrix r) { // 矩陣相乘
int dim = q[0].size();
if (dim != r.size()) return F_matrix(0);
int m = q.size();
int n = r[0].size();
F_matrix x(m, F_vector(n));
for (int i = 0; i < m; i ++) {
for (int j = 0; j < n; j ++) {
double sum = 0;
for (int k = 0; k < dim; k++) {
sum += q[i][k] * r[k][j];
}
x[i][j] = sum;
}
}
return x;
}
void jacobiSVD(F_matrix &A) {// A = USVᵗ, A -> U[σ]V
int mm = A.size();
int nn = A[0].size();
F_matrix aU = A;// 初始化 L 左矩陣
F_matrix aV = eye(nn);// 初始化 V 正交矩陣
if (mm < nn) {
aU.resize(nn);
for (int i = mm; i < nn; i++) { aU[i].resize(nn); }
mm = nn;
}
int iterations = 30;
while (iterations -- > 0) {
double converge = 0;
for (int i = 1; i < nn; i ++) {
for (int j = 0; j < i; j ++) {
double a = 0;
double b = 0;
double c = 0;
for (int k = 0; k < mm; k ++) {
a += aU[k][j] * aU[k][j];// a = L̂ⱼ•L̂ⱼ
b += aU[k][i] * aU[k][i];// b = L̂ᵢ•L̂ᵢ
c += aU[k][j] * aU[k][i];// c = L̂ⱼ•L̂ᵢ = L̂ᵢ•L̂ⱼ
}
if ( a < tol || b < tol) { continue; }
double s = (c < 0) ? -c / sqrt(a*b) : c / sqrt(a*b);
if (s < tol) { continue; }
if (s > converge) converge = s;
F_matrix J = eye(nn);// 初始化旋轉矩陣 Jⱼᵢ(θ) = I
s = (b - a) / (2 * c);// 只要算出 cosθ, sinθ, 不用算出旋轉角 θ
double t = (s < 0) ? -1.0/(-s + sqrt(1 + s*s)) : 1.0/(s + sqrt(1 + s*s));
c = 1 / sqrt(1 + t * t);
s = t * c;
J[j][j] = c;
J[i][j] = -s;
J[j][i] = s;
J[i][i] = c;
aU = product(aU, J);// 更新 L 左矩陣
aV = product(aV, J);// 更新 V 正交矩陣
}
}
if (converge < 1e-9) break; // todo: break when stable
}
printf("Singular value:\nσ:\t");
for (int i = 0; i < nn; i ++) {
double sum = 0;
for (int k = 0; k < mm; k ++) {
sum += aU[k][i] * aU[k][i]; // Σ(L̂ᵢ[k])²
}
double sigma = sqrt(sum);// Singular value σᵢ = ||L̂ᵢ||₂ = √[Σ(L̂ᵢ[k])²]
if (sigma > 1e-9) {// sigma != 0
for (int k = 0; k < mm; k ++) {// ûᵢ = L̂ᵢ/σᵢ, that is why σ name "singular value"
aU[k][i] /= sigma;
}
}
printf("%-10.4f\t", sigma);
}// todo： 未完成排序, 當 σ = 0 需刪除, 對應的 û/v̂ 向量也要跟著移除.
printf("\n");
matrixShow(aU);
matrixShow(aV);
}
F_matrix arrayclone(double *array, int w, int h) {
F_matrix temp(h, F_vector(w));
for (int i = 0; i < h; i++) {
for (int j = 0; j < w; j++) { temp[i][j]= *array ++; }
}
return temp;
}
int main(int argc, char **argv) {
double a[] = {
1, 2, 3, 4, 5,
6, 7, 8, 9,10
};
F_matrix B = arrayclone(a, 5, 2);
F_matrix C = transpose(B);
matrixShow(B);
jacobiSVD(B);
printf("\n======================\n");
matrixShow(C);
jacobiSVD(C);
}
