如何用C++实现SVD算法

SVD（奇异值分解）是一种常用的矩阵分解方法，可以用于数据降维、压缩、去噪等方面。以下是用C++实现SVD算法的基本步骤：

读入需要分解的矩阵A，假设其大小为m x n。
计算A的转置矩阵AT，大小为n x m。
计算A和AT的乘积，得到n x n的对称矩阵ATA。
对ATA进行特征值分解，得到特征值和特征向量。
将特征向量按特征值从大到小排序，得到n x n的正交矩阵U。
计算矩阵V=ATU，得到n x n的正交矩阵V。
计算矩阵S，大小为m x n，其中S的对角线元素为ATA的特征值的平方根，其余元素为0。
将S的大小调整为m x n，同时将U和V的大小调整为m x m和n x n。
输出U、S和V，即为矩阵A的SVD分解结果。

以下是C++代码实现SVD算法的基本框架：

#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>

using namespace std;

// 定义矩阵类型
typedef vector<vector<double>> Matrix;

// 计算矩阵转置
Matrix transpose(const Matrix &A) {
    int m = A.size(), n = A[0].size();
    Matrix AT(n, vector<double>(m));
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            AT[j][i] = A[i][j];
        }
    }
    return AT;
}

// 计算矩阵乘积
Matrix multiply(const Matrix &A, const Matrix &B) {
    int m = A.size(), n = A[0].size(), p = B[0].size();
    Matrix C(m, vector<double>(p));
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < p; j++) {
            for (int k = 0; k < n; k++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
    return C;
}

// 计算矩阵的Frobenius范数
double norm(const Matrix &A) {
    double sum = 0;
    for (auto row : A) {
        for (auto x : row) {
            sum += x * x;
        }
    }
    return sqrt(sum);
}

// 计算矩阵的特征值和特征向量
void eigen(const Matrix &A, vector<double> &eigval, Matrix &eigvec) {
    // TODO: 实现特征值分解算法
}

// 计算SVD分解
void svd(const Matrix &A, Matrix &U, Matrix &S, Matrix &V) {
    int m = A.size(), n = A[0].size();
    Matrix AT = transpose(A);
    Matrix ATA = multiply(AT, A);

    // 计算ATA的特征值和特征向量
    vector<double> eigval;
    Matrix eigvec;
    eigen(ATA, eigval, eigvec);

    // 将特征向量按特征值从大到小排序
    vector<pair<double, vector<double>>> pairs;
    for (int i = 0; i < n; i++) {
        pairs.push_back({eigval[i], eigvec[i]});
    }
    sort(pairs.rbegin(), pairs.rend());

    // 构造U和V矩阵
    U.resize(m, vector<double>(m));
    V.resize(n, vector<double>(n));
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < m; j++) {
            U[i][j] = pairs[j].second[i];
        }
    }
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            V[i][j] = pairs[j].second[i];
        }
    }

    // 计算S矩阵
    S.resize(m, vector<double>(n));
    for (int i = 0; i < n; i++) {
        double eig = sqrt(pairs[i].first);
        S[i][i] = eig;
    }

    // 调整S、U和V的大小
    if (m > n) {
        S.resize(n, vector<double>(n));
        U.resize(m, vector<double>(n));
    } else if (m < n) {
        S.resize(m, vector<double>(m));
        V.resize(n, vector<double>(m));
    }
}

int main() {
    // 读入矩阵A
    int m, n;
    cin >> m >> n;
    Matrix A(m, vector<double>(n));
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            cin >> A[i][j];
        }
    }

    // 计算SVD分解
    Matrix U, S, V;
    svd(A, U, S, V);

    // 输出结果
    cout << "U:" << endl;
    for (auto row : U) {
        for (auto x : row) {
            cout << x << " ";
        }
        cout << endl;
    }
    cout << "S:" << endl;
    for (auto row : S) {
        for (auto x : row) {
            cout << x << " ";
        }
        cout << endl;
    }
    cout << "V:" << endl;
    for (auto row : V) {
        for (auto x : row) {
            cout << x << " ";
        }
        cout << endl;
    }

    return 0;
}