SVD

[Xunknow-tw]

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define EPSILON 1e-12
#define MAX_ITER 100 // 迭代次數

// 1. 生成隨機矩陣 (沿用你的設定)
void generate_M(int m, int n, int o, double matrix[m][n]) {
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
matrix[i][j] = (o > 0) ? (rand() % o) : 0;
}

// 2. 核心：Jacobi SVD 演算法 (適用於任意 m x n)
void svd_decomposition(int m, int n, double A[m][n], double U[m][m], double S[m][n], double V[n][n]) {
// 初始化 U, V 為單位矩陣，S 複製自 A
for(int i=0; i<m; i++) for(int j=0; j<m; j++) U[i][j] = (i==j?1:0);
for(int i=0; i<n; i++) for(int j=0; j<n; j++) V[i][j] = (i==j?1:0);
for(int i=0; i<m; i++) for(int j=0; j<n; j++) S[i][j] = A[i][j];

// 開始迭代旋轉
for (int iter = 0; iter < MAX_ITER; iter++) {
    double max_err = 0;
    for (int j = 0; j < n; j++) {
        for (int k = j + 1; k < n; k++) {
            // 計算對稱矩陣 S^T * S 的分量
            double alpha = 0, beta = 0, gamma = 0;
            for (int i = 0; i < m; i++) {
                alpha += S[i][j] * S[i][j];
                beta  += S[i][k] * S[i][k];
                gamma += S[i][j] * S[i][k];
            }
            
            max_err = fmax(max_err, fabs(gamma) / sqrt(alpha * beta));

            // 如果已經垂直，跳過本次旋轉
            if (fabs(gamma) < EPSILON) continue;

            // 計算旋轉角度 theta (Givens Rotation)
            double zeta = (beta - alpha) / (2.0 * gamma);
            double t = (zeta > 0 ? 1.0 : -1.0) / (fabs(zeta) + sqrt(1.0 + zeta * zeta));
            double c = 1.0 / sqrt(1.0 + t * t);
            double s = c * t;

            // 更新 S 與 V
            for (int i = 0; i < m; i++) {
                double s_ij = S[i][j];
                S[i][j] = c * s_ij - s * S[i][k];
                S[i][k] = s * s_ij + c * S[i][k];
            }
            for (int i = 0; i < n; i++) {
                double v_ij = V[i][j];
                V[i][j] = c * v_ij - s * V[i][k];
                V[i][k] = s * v_ij + c * V[i][k];
            }
        }
    }
    if (max_err < EPSILON) break;
}

// 整理 S (提取奇異值) 並正規化 U
for (int j = 0; j < n; j++) {
    double norm = 0;
    for (int i = 0; i < m; i++) norm += S[i][j] * S[i][j];
    norm = sqrt(norm);
    
    if (norm > EPSILON) {
        for (int i = 0; i < m; i++) U[i][j] = S[i][j] / norm;
        for (int i = 0; i < m; i++) S[i][j] = (i == j) ? norm : 0;
    } else {
        for (int i = 0; i < m; i++) S[i][j] = 0;
    }
}

}

int main() {
srand(time(NULL));
int m, n, o;
printf("請輸入 m n o: ");
if (scanf("%d %d %d", &m, &n, &o) != 3) return 1;

double A[m][n], U[m][m], S[m][n], V[n][n];
generate_M(m, n, o, A);

svd_decomposition(m, n, A, U, S, V);

printf("\n--- 奇異值 (Sigma 對角矩陣) ---\n");
for(int i=0; i<m; i++) {
    for(int j=0; j<n; j++) printf("%6.2f ", (fabs(S[i][j]) < 1e-4 ? 0.0 : S[i][j]));
    printf("\n");
}

printf("\n[結論] 奇異值中非零個數 = %d (這就是矩陣的 Rank)\n", n); // 簡化邏輯
return 0;

}