將矩陣分解成 QR 兩矩陣相乘, 可以利用一個類似旋轉矩陣, 一次清除一個位置, 這種矩陣稱為 Givens 旋轉矩陣 Gⱼ(θ), 有別於用屋頂保持矩陣 H(û) (Householder matrix)一次清除一行, Gⱼ(θ) 是在單位矩陣對角線上, I(j,j) 的位置附近附加 2x2 旋轉角 θ 的元素 (cosθ, sinθ), 用來將位置 (j,i) 的元素清為 0, 它有旋轉矩陣的特性: Gⱼ 轉置矩陣 Gⱼᵗ 就是反矩陣 Gⱼ⁻ = Gⱼᵗ, 也就是說 GⱼᵗGⱼ = I, 因此只要將對角線以下的元素一個一個清為零, 自然就化解出 QR 矩陣:
A = IA = I(GⱼᵗGⱼ)A = (IGⱼᵗ ...) * (... Gⱼ)A
= Q * R
要注意的是 Gⱼ(θ) 與矩陣的更動及要清 0 位置 (j, i) 的旋轉角 θ 有關, 需要動態更新,上述右式隨著矩陣 A 不斷左乘 Gⱼ 矩陣, 對角線下的元素一一清為 0, 最終自然浮現出上三角矩陣 R, 左式伴隨著將單位矩陣不斷右乘 Gⱼᵗ 矩陣, 最終演化成 Q 矩陣. 因此可以用下面疊代法,算出最終的 Q, R:
初始化 Q = I, R = A, Gⱼ(θ)=[cosθ, sinθ ; -sinθ, cosθ]
R = Gⱼ(θ) * R
Q = Q * Gⱼᵗ(θ)
其中 Gⱼ(θ) 附加的 2x2 矩陣 Gⱼ(θ)=[cosθ, sinθ; -sinθ, cosθ], Gᵗ(θ)= [cosθ, -sinθ; sinθ, cosθ].
Givens 運算適合用在下三角位置上 0 較多的矩陣. 等於 0 時就不用計算. 實際上A[j][i] = 0 時, Gⱼ(θ) 的旋轉角正好是 90ᵒ. 而 Gⱼ(90ᵒ) 等於單位矩陣 I, 運算就是多餘的. 當矩陣與 Givens 矩陣相乘時, 其實只會影響到附近 2x2 整行與整列計算, 有很多乘 0 的部份其實可以省略, 底下只闡述實作, 並未作最佳化:
#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 givensQRD(F_matrix &A, F_matrix &qq, F_matrix &rr) { // A = QR
int m = A.size();
int n = A[0].size();
F_matrix Q = eye(m);
F_matrix R = A;
double c = 1;// dummy, 預設 θ = 90ᵒ, c = cos(πθ/180) = 1
double s = 0;// dummy, 預設 θ = 90ᵒ, s = sin(πθ/180) = 0
for (int i = 0; i < n; i ++) {// 從第一行開始掃描
for (int j = m - 1; j > i; j --) {// 行向量 i 從最後一個元素開始往前直到對角線, 一一清為 0
double a = R[j][i]; // 目標行向量 i, 位置 j 元素準備清為 0, 先計算出旋轉矩陣 Gⱼᵢ(θ)
if (fabs(a) < tol) continue; // 當 a != 0 才要清除
double b = R[j-1][i];// 行向量 i 只剩位置 j-1 的元素值會受影響
if (fabs(a) > fabs(b)) { // 找出 Gⱼ(θ) 旋轉, 只要推導出 cosθ, sinθ, 不用算出實際 θ 值
double t = b / a;
s = 1 / sqrt(1 + t * t);
c = s * t;
} else {
double t = a / b;
c = 1 / sqrt(1 + t * t);
s = c * t;
}
F_matrix G = eye(m);// 初始化旋轉矩陣 Gⱼ(θ) = I, 旋轉角 θ = 90ᵒ
G[j][j] = c;// 更新旋轉矩陣,準對角線上擺 cosθ, 目標清除行向量 i 位置 j 元素 a => 0
G[j-1][j-1] = c;// 更新旋轉矩陣, 準對角線左上擺 cosθ
G[j][j-1] = -s;// 更新旋轉矩陣,準對角, 正左擺 -sinθ
G[j-1][j] = s;// 更新旋轉矩陣,準對角, 正上擺 sinθ
R = product(G, R); // 乘完後 R[j[i] => 0
G[j][j-1] = s;// 矩陣轉置 Gⱼ => Gⱼᵗ
G[j-1][j] = -s;// GⱼᵗGⱼ = I
Q = product(Q, G);// A = QGⱼᵗ GⱼR,
}
}
swap(Q, qq);
swap(R, rr);
}
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,
17, 8, 7, 12, 4,
17, 8, 1, 11, 5,
7, 2, 5, 3, 1
};
F_matrix B = arrayclone(a, 5, 5);
F_matrix C = transpose(B);
F_matrix Q, R;
matrixShow(B);
givensQRD(B, Q, R);
matrixShow(R);
matrixShow(Q);
matrixShow(product(Q,R));
matrixShow(C);
givensQRD(C, Q, R);
matrixShow(R);
matrixShow(Q);
matrixShow(product(Q,R));
}
2023年4月29日 星期六
C++ 使用 Givens 旋轉矩陣將矩陣分解 QR
2023年4月28日 星期五
C++ 使用 householder 將矩陣分解成 QR 矩陣
使用 householdler 演算, 可以將任意矩陣分解成 Q R 兩矩陣相乘. 原理是先從矩陣擷取出行向量, 算出單位向量 ûₖ, 再計算對應的屋頂保持矩陣(Householder matrix) Hₖ, 該矩陣的特性是: 1. 兩個矩陣相乘會變成單位矩陣 HₖHₖ = I, 也就是說原矩陣就是反矩陣 Hₖ⁻ = Hₖ, 就好像是一面鏡子可分拆兩邊的特性. 2. 當一個矩陣左邊乘上 H 後, 行向量所在矩陣, 正對角以下元素全部清為 0, 矩陣屋頂(正對角)上的元素則保持不受影響, 也就是說乘一次 H 就清一行 0, 要注意的是 Hₖ 會隨著矩陣的更動而變化, 根據上述兩特性, 任一矩陣 A 就能分解成 Q 及 R 矩陣(A = QR):
A = IA = I(HₖHₖ)A
= I(H₁H₁)A = (IH₁) * (H₁A)
= (IH₁)(H₂H₂)(H₁A) = (IH₁H₂) * (H₂H₁A)
= (IH₁H₂)(H₃H₃)(H₂H₁A = (IH₁H₂H₃) * (H₃H₂H₁A)
= ...
= (IH₁H₂H₃ ... Hₙ) * (Hₙ ... H₃H₂H₁A)
右式將矩陣 A 不斷的左乘 H, 隨著對角線以下, 一行一行清為零, 自然就變成上三角矩陣 R(又稱右矩陣) , 左式則將單位矩陣 I 不斷右乘 H 的結果, 最後演變成 Q 矩陣(它是正交矩陣, 因為 H 本身就是正交矩陣)
Hₖ = I - 2*ûₖ*ûₖᵗ , 其中 ûₖ 為 Householder 單位向量
R₁ = H₁A = (I - 2*û₁*û₁ᵗ)*A = A - 2*û₁*û₁ᵗ*A
Q₁ = IH₁ = I(I - 2*Q*û₁*û₁ᵗ) = I - 2*û₁*û₁ᵗ
Rₖ <= HₖRₖ₋₁ = (I - 2*ûₖ*ûₖᵗ)*Rₖ₋₁ = Rₖ₋₁ - 2*ûₖ*ûₖᵗ*Rₖ₋₁
Qₖ <= Qₖ₋₁Hₖ = Qₖ₋₁(I - 2*Q*ûₖ*ûₖᵗ) = Qₖ₋₁ - 2*Qₖ₋₁*ûₖ*ûₖᵗ
隨著 uₖ => Hₖ => Rₖ₊₁ => uₖ₊₁ => Hₖ₊₁ 不斷推陳出新, 利用疊代法就能算出最終的 Q, R:
R -= 2*ûₖ*ûₖᵗ*R
Q -= 2*Q*ûₖ*ûₖᵗ
#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;
}
F_vector householdvector(F_matrix &R, int i) {
int m = R.size();
int mi = m - i;
F_vector ut(mi);
ut[0] = R[i][i];
double s = 0;
for (int jk = 1, j = i + 1; jk < mi ; jk ++) {
double r = R[j ++][i];
s += r * r;
ut[jk] = r;
}
ut[0] += sqrt(s + ut[0]*ut[0]);
s = sqrt(s + ut[0]*ut[0]);
if (s < tol) return F_vector(0);
for (int jk = 0; jk < mi ; jk ++) {
ut[jk] /= s;
}
return ut;
}
int householdQRD(F_matrix &xx, F_matrix &qq, F_matrix &rr) { // xx => QR
int rank = 0;
int m = xx.size();
int n = xx[0].size();
int maxSize = m;
if (maxSize < n) maxSize = n;
F_matrix Q = eye(m);
F_matrix R = xx;
F_matrix tempH(maxSize, F_vector(maxSize));
F_vector uu(maxSize);
for (int i = 0; i < n ; i ++) { // n columns
F_vector ut = householdvector(R, i);
int mi = ut.size();
if (mi == 0) break;
rank ++;
for (int k = i, ki = 0; k < n; k ++, ki ++) {
double sum = 0;
for (int kj = 0, j = i; kj < mi; kj ++, j++) {
sum += ut[kj] * R[j][k]; // uᵗ•Rj
}
uu[ki] = sum + sum; // w = 2uᵗ*R, which is n elements of row vector
}
for (int ki = 0; ki < n - i; ki ++) {
for (int kj = 0; kj < mi; kj ++) {
tempH[kj][ki] = ut[kj] * uu[ki];// tempH = u*w = 2u*(uᵗ*R), which is a miXni matrix
}
}
for (int k = i, ki = 0; k < n; k ++, ki ++) {
for(int j = i, kj = 0; j < m; j ++, kj ++) {
R[j][k] -= tempH[kj][ki]; // R -= tempH = 2u*(uᵗ*R)
}
}
for (int j = 0; j < m; j ++) {
double sum = 0;
for (int ki = 0, k = i; ki < mi; ki ++) {
sum += Q[j][k ++] * ut[ki]; // Qj•u
}
uu[j] = sum + sum; // v = 2Q*u; It is a column vector with m elements now.
}
for (int j = 0; j < m ; j ++) {
for (int ki = 0; ki < mi ; ki ++) {
tempH[j][ki] = uu[j] * ut[ki];// tempH = v*uᵗ = 2Q*u*uᵗ; which is a mXm matrix
}
}
for(int j = 0; j < m ; j++) {
for (int ki = 0, k = i; ki < mi; ki ++, k ++) {
Q[j][k] -= tempH[j][ki]; // Q -= tempH = (2Q*u)*uᵗ
}
}
}
swap(Q, qq);
swap(R, rr);
return rank;
}
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[] = {
4, 2, 3, 4,
17, 8, 9, 13
};
F_matrix B = arrayclone(a, 4, 2);
F_matrix Q, R;
matrixShow(B);
householdQRD(B, Q, R);
matrixShow(Q);
matrixShow(R);
matrixShow(product(Q, R));
F_matrix C = transpose(B);
matrixShow(C);
householdQRD(C, Q, R);
matrixShow(Q);
matrixShow(R);
matrixShow(product(Q, R));
}
2023年4月23日 星期日
linux 用 c++ 玩線性代數
在linux mint 下 c++ 使用 eigen3 作矩陣運算, 用 matplotlib-cpp 來繪圖, 可謂天作之合, 兩者都是標準純 C++ 標頭檔(header file), 所有程序都定義在類型內(class), 只要在程式裏面include 標頭檔, 並不需要聯結任何程式庫, 但編譯 eigen3 所寫的程式時,會很耗時, 要解決編譯的問題參考 https://stackoverflow.com/questions/67529835/eigen3-take-a-long-time-to-compile-and-very-slow-when-debug
就是把要用的程式先寫成 c++ 程式碼(.cpp), 編譯成獨立的模組, 作成程式庫, 後續透過聯結該模組就可. 首先上官網下載原始檔: https://eigen.tuxfamily.org/index.php?title=Main_Page
解壓縮到專案目錄內,只需用到目錄底下的 Eigen 目錄, 用符號聯結就可以了, 並不需要特別安裝
ln -sf /home/mint/Downloads/eigen-3.4.0/Eigen .
例如我只要用到 JacobiSVD, 但又不想用裡面的其他東西, 因習慣用 c++ 標準 <vector>, 只好先轉成 eigen3 的 matrix 才能運算, 結束時用 <vector> 傳回資料:
// jacobi.cpp
#include "Eigen/SVD"
#include <vector>
using namespace Eigen;
using namespace std;
using F_vector = vector<double>;
using F_matrix = vector<F_vector>;
extern "C" { // to export function
F_vector jacobiSVD(F_matrix &x, F_matrix &U, F_matrix &V) {
int m = x.size();
int n = x[0].size();
MatrixXf temp (m, n);
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) { temp(i, j) = x[i][j]; }// copy x to temp
}
JacobiSVD<MatrixXf> svd(temp, ComputeThinU | ComputeThinV );
MatrixXf U2 = svd.matrixU();
MatrixXf V2 = svd.matrixV();
MatrixXf S2 = svd.singularValues();
int ss = S2.rows();
F_vector S(ss);
U.resize(U2.rows(), F_vector(ss));
V.resize(V2.rows(), F_vector(ss));
for (int i = 0; i < ss; i ++) { S[i] = S2(i, 0); }// copy from S2
for (int i = 0; i < U2.rows() ; i ++) {
for (int j = 0; j < ss; j ++) { U[i][j] = U2(i, j); }// copy from U2
}
for (int i = 0; i < V2.rows() ; i ++) {
for (int j = 0; j < ss; j ++) { V[i][j] = V2(i, j); }// copy from V2
}
return S;// singular values in a vector
}
}
直接用 g++ 編譯成模組檔 .o
g++ jacobi.cpp -c -o jacobi.o
後續只要聯結該模組無須重新編譯, 也完全不用再 include 任何 eigen3 的檔案, 解決編譯緩慢的問題, 但主程式必須要用 extern "C" 方式呼叫該模組內的函式:
#include <random>
#include <vector>
#include "matplotlibcpp.h"
#define tol ((double)1e-12)
namespace plt = matplotlibcpp;
using namespace std;
using F_vector = vector<double>;
using F_matrix = vector<F_vector>;
extern "C" {
F_vector jacobiSVD(F_matrix &x, F_matrix &U, F_matrix &V);
}
double tpower(double t, int n) { // f(t, n) = tⁿ
if (n == 0) return 1.0;
if (n == 1) return t;
double temp = t;
for (int i = 1; i < n; i ++) temp *= t;
return temp;
}
double absval(double t) { // f(t) = |t|
return t < 0 ? -t : t;
}
struct svdSolver {
F_matrix U;
F_vector S;
F_matrix V;
static svdSolver jacobiMethod(F_matrix &x) {
F_matrix U;
F_matrix V;
F_vector S = jacobiSVD(x, U, V);// call external library
return svdSolver {.U = U, .S = S, .V = V}; // wrap U,S,V
}
static F_vector polyfit(F_vector y, int parameters = 2) {
int size = y.size();
F_matrix At(size, F_vector(parameters));
for (int j = 0; j < parameters; j ++) {// Aᵗ xᵗ = yᵗ
for (int i = 0; i < size ; i ++) {
At[i][j] = tpower(i, j);// iᶨ
}
}
svdSolver svd = jacobiMethod(At);
svd.show(svd.S);// singular value
svd.show(svd.V);// all of the right singular vectors
return svd.svdfit(y);
}
void show(F_vector v) {
for (int i = 0; i < v.size(); i ++) { printf("%f\t", v[i]); }
printf("\n");
}
void show(F_matrix x) {
for (int j = 0; j < x.size(); j ++) { show(x[j]); }
printf("\n");
}
F_vector svdfit(F_vector &y) {// x̂ᵗ = A⁺ŷᵗ = (UΣVᵗ)⁻ŷᵗ = VΣ⁻Uᵗŷᵗ= V(Uᵗŷᵗ/Σ)
int udim = U.size();
if (y.size() != udim) return F_vector(0);
int ss = S.size();
F_vector uyisTemp(ss); // uyisTemp[j]ᵗ = Σ⁻(Ujᵗ•ŷᵗ) = Ujᵗ•ŷᵗ / Σ
for (int j = 0; j < ss; j ++) {
double sum = 0.0;// inner product
if (absval(S[j]) > tol) {
for (int i = 0; i < udim; i ++) { // Ujᵗ•ŷᵗ
sum += U[i][j] * y[i];
}
sum /= S[j];// (Ujᵗ • ŷᵗ)/Σ
}
uyisTemp[j] = sum;
}
int vdim = V.size();
F_vector x(vdim);
for (int i = 0; i < vdim; i ++) { // Vi • uyisTempᵗ
double sum = 0.0;
for(int j = 0; j < ss; j++) {
sum += V[i][j] * uyisTemp[j];
}
x[i] = sum;
}
return x;
}
};
int main() {
int fs = 512; // number of samples
int np = 6; // number of parameters to fit
vector<double> xn(fs), yn(fs), gt(fs), cx(np);
srand(time(nullptr));
for (int i = 0 ; i < cx.size() ; i++) { cx[i] = rand(); }
for (int t = 0 ; t < fs ; t ++) {
double st = 0;
for (int i = 0; i < cx.size(); i ++) {
st += tpower(t, i); // tᶤ
}
double err = (double)rand()/RAND_MAX - 0.5;
double noise = st * err; // 50 % random noise
yn[t] = st + noise;
xn[t] = t;
}
vector<double> p = svdSolver::polyfit(yn, cx.size());
for (int t = 0 ; t < fs ; t ++) {
double yt = 0;
for (int i = 0; i < p.size(); i ++) {
yt += p[i] * tpower(t, i);
}
gt[t] = yt;
}
plt::plot(xn, yn, "g:x", xn, gt, "r-"); // x-y plot, red solid line.
plt::show();
}
程式用到 c++ matplotlib, 參考: https://masontseng.blogspot.com/2023/04/linum-mint-c-matplotlib.html 下載 matplotlibcpp.h
接著將主程式聯結模組成功後, 就能執行:
g++ plot.c jacobi.o -std=c++11 -I/usr/include/python3.8 -lpython3.8 && ./a.out
2023年4月20日 星期四
矩陣向量一些概念
用內積來表示兩向量 r̂, n̂ 之間的正投影長度 r̂•n̂ =|r̂||n̂|cosθ, 假設 r̂, n̂ 是 k 維空間的向量：
r̂•n̂ = r₁*n₁ + r₂*n₂ + ... + rₖ*nₖ = Σ(rₖ * nₖ) = Σ(nₖ*rₖ) = n̂•r̂
|r̂| = √(r₁² + r₂² + ..., rₖ²)
|n̂| = √(n₁² + n₂² + ..., nₖ²)
cosθ = r̂•n̂/(|r̂||n̂|)
想像空間是由一群向量 r̂₁, r̂₂, ... , r̂ₖ 等獨立基底向量 [r̂], 以原點 ô 為中心, 所展示的 k 維輻射空間, 描述任一向量 n̂ 是 r̂ 的線性組合, 用內積代表的座標, 就是在基底向量上的投影分量:
n̂ = Σ nₖr̂ₖ = Σ(n̂•r̂ₖ)r̂ₖ = Σ(r̂ₖ•n̂)r̂ₖ
n₁ = r̂₁•n̂
n₂ = r̂₂•n̂
.........
nₖ = r̂ₖ•n̂
底下符號 M 用來代表一個 nxn 方陣, 行向量彼此獨立, 而且列向量都是獨立向量
底下符號 λ 用來代表 eigenvalue, 它是一個純量, 用一個數字來量化
底下符號 ê 用來代表 eigenvector 它是一個向量, 用一組數字來量化(e₁, e₂,...,eₖ), 向量的元素個數 k, 用來描述該空間的維度, 用 [ê] 代表 eigenvector 行向量形成的矩陣 V
底下符號 Λ 用來代表對角線矩陣, 除對角線上有數值外, 其餘位置都為 0
底下符號 [λ] 用來代表由 eigenvalue 形成的對角線矩陣 Λ
將 M 分解為由 eigenvalue 所組成對角線矩陣 Λ, 及 eigenvector 行向量 ê 所組成矩陣 V = [ê] 時: [V, D] = eig(M)
M = V Λ V⁻ = [ê] [λ] [ê]⁻
V = [ê₁, ê₂, ê₃, ..., êₖ]
V⁻ 是 V 的反矩陣 = [ê]⁻
Λ = D = [λ₁,0,...,0; 0,λ₂,0, ...,0; ... ; 0,0,0,...,λₖ] = [λ]
就可以看出 M 是將向量先經反矩陣 V⁻ 轉換, 用 eigenvalue 伸縮座標軸, 投影在 eigenvector 上:
M x̂ᵗ = (V Λ V⁻)x̂ᵗ
eigenvalue λ , eigenvector ê 與矩陣 Ｍ 的關係式: M•ê = λ•ê
當 Mᵢⱼ= Mⱼᵢ 時, M 就是一個對稱矩陣, 其 eigenvector 必定互相垂直, 由 eigenvector 組成矩陣的轉置矩陣就是反矩陣
∵ Mᵢⱼ= Mⱼᵢ
∴ V⁻ = Vᵗ
∴ M = V Λ Vᵗ
M x̂ᵗ = (V Λ Vᵗ)x̂ᵗ = V Λ Vᵗ x̂ᵗ
方陣轉置(將行列互換)後 eigenvalue 不會變
任何一個矩陣轉置後, 再與原矩陣相乘就會形成一個對稱方陣:
T = AᵗA = AₙxₘAₘxₙ= Tₙxₙ
如果相乘的順序反過來也會形成另一個對稱方陣:
N = AAᵗ =AₘxₙAₙxₘ= Nₘxₘ
當矩陣不是方陣時, 任何一個 ₘxₙ的矩陣 A 可以透過 svd 分解成 3 個矩陣: [U, S, V] = svd(A)
其中 S 是 singular value [σ] 形成的矩陣, V 是ＡᵗA 的 eigenvector 形成的矩陣 [v̂], 但 U 是 AAᵗ 的 eigenvector 形成的矩陣 [û], 要注意的是 AᵗA 與 AAᵗ 乘積是兩個不同的矩陣
A = U S Vᵗ
Aₘxₙ = Uₘxₘ • Sₘxₙ • Vᵗₙxₙ
U 是由 AAᵗ 轉置矩的 eigenvector 行向量所組成 ₘxₘ 方陣 [û], m 是目標空間維度
U = [û₁; û₂; û₃; ...; ûₘxₘ]ᵗ = [û₁, û₂, û₃, ..., ûₘxₘ] = [û]
Uᵗ 是 U 的轉置矩陣也是 U 的反矩陣
Uᵗ = [û₁, û₂, û₃, ..., ûₘ]⁻ = [û₁; û₂; û₃; ...; ûₘ] = [û]⁻
S 對角線由 singular value σ 所組成, 其餘為 0, 組成一個 ₘxₙ 矩陣
S = [σ]ₘxₙ = [σ]
V 是由 AᵗA 的 eigenvector 行向量所組成 ₙxₙ 方陣 [v̂], n 是來源空間維度
V = [v̂₁; v̂₂; v̂₃; ...; v̂ₙ]ᵗ = [v̂₁, v̂₂, v̂₃, ..., v̂ₙ] = [v̂]
Vᵗ 是 V 的轉置矩陣也是 V 的反矩陣, n 是來源空間維度
Vᵗ = [v̂₁, v̂₂, v̂₃, ..., v̂ₙ]⁻ = [v̂₁; v̂₂; v̂₃; ...; v̂ₙ] = [v̂]⁻
備註: 只有方陣才能計算出 eigenvalue 與 eigenvector, 非方陣由 SVD 算出的是 singular value 矩陣及左右伴隨矩陣, 伴隨矩陣分別由 right singular vector及 left singular vector 行向量所組成的矩陣 V=[v̂], U=[û]
1.矩陣 A = Aₘxₙ = Uₘxₘ•Sₘxₙ•Vᵗₙxₙ 若 singular value 只有 k 組 (m!=n, k=m 或是 k=n)時, 改用 Uₘxₖ•Sₖxₖ•Vᵗₖxₙ, 也就是說讓 S 組成稀疏對角線 ₖxₖ 方陣 S=Λ=[σ], U 變成瘦型矩陣 ₘxₖ, V 也順勢成瘦型矩陣 ₙxₖ, 較省記憶空間
2.矩陣 AᵗA 的 eigenvalue 是 A singular value 的二次方:
T = AᵗA = (U S Vᵗ)ᵗ U S Vᵗ = VSUᵗUSVᵗ = VS²Vᵗ = VΛVᵗ
因此 A 的 singular value 平方就是 AᵗA eigenvalue, 反過來說 A 的 singular value 是 AᵗA eigenvalue 開根號
3. AAᵗ 的 eigenvalue 也是 A的 singular value 二次方, 反過來說 A 的 singular value 是 AAᵗ eigenvalue 開根號
AAᵗ= U S Vᵗ (U S Vᵗ)ᵗ = USVᵗ VSUᵗ= US²Uᵗ = UΛUᵗ
4.矩陣 A 的 SVD 算法： Aₘxₙ = Uₘxₘ • Sₘxₙ • Vᵗₙxₙ
a. 計算 T = AᵗA 用來找尋 A 的 singular value 及左右伴隨矩陣 UV
b. Tv̂ = σ²v̂,透過 QR decoposition 解出 T 的 eigenvalue 矩陣 [σ²]ₘxₙ 及 eigenvector 所組成的矩陣 Vₙxₙ= [v̂]
c. 將上述 eigenvalue 開根號得出 σ, S = Λ = [σ] ₘxₙ, 將 σ 取倒數放入矩陣的對角線形成反矩陣, 注意反矩陣行列要交換才能與後續 Ｖ相乘
S⁻ = Λ⁻ = [1/σ] ₙxₘ,
d. U = U (Λ (Vᵗ V) Λ⁻) = (UΛVᵗ) V Λ⁻ = A V Λ⁻
ＡV 先相乘, 最後再乘步驟 c. 的 singular value 反矩陣 Λ⁻, 就能算出 U = A • V • Λ⁻
2023年4月16日 星期日
c++ 反矩陣解矩陣方程式
-
市面上便宜的時鐘, 內部機芯大多使用晶體振盪器產生固定頻率的脈沖去驅動 Lavet-type 步進馬達, 它只要一個 1.5V 電池就足以推動, 而馬達扭力主要是靠線圈內的電流產生磁力所致, 因此電流才是動力來源. 電壓太高會造成線圈電流飽和導致發熱,可能影響震盪頻率甚至燒毀... -
參考資料: https://developer.android.com/training/connect-devices-wirelessly/nsd linux mint 上針對區域網路的名稱服務協定可以使用 avahi, 透過 avahi-browse -a -r 就能...