矩阵SVD分解计算

对于有数学库的时候，进行矩阵相关计算还是不复杂，但是没有数学库就很麻烦，利用算法实现了矩阵奇异值分解。

void decompose(const std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& U, std::vector<double>& S,
        std::vector<std::vector<double>>& V) 
{
        if (A.empty() || A[0].empty()) 
        {
            throw std::invalid_argument("Matrix A cannot be empty");
        }

        size_t m = A.size();
        size_t n = A[0].size();

        if (m < n) 
        {
            throw std::invalid_argument("Number of rows must be >= number of columns in SVD");
        }

        U = A;
        S.resize(n);
        V.assign(n, std::vector<double>(n, 0.0));
        
        std::vector<double> rv1(n);
        double g = 0.0, scale = 0.0, anorm = 0.0;

        for (size_t i = 0; i < n; ++i) {
            size_t l = i + 1;
            rv1[i] = scale * g;
            g = scale = 0.0;

            if (i < m) {
                for (size_t k = i; k < m; ++k) scale += std::abs(U[k][i]);
                if (scale != 0.0) {
                    for (size_t k = i; k < m; ++k) {
                        U[k][i] /= scale;
                        g += U[k][i] * U[k][i];
                    }
                    double f = U[i][i];
                    g = -sign(std::sqrt(g), f);
                    double h = f * g - g * g;
                    U[i][i] = f - g;
                    for (size_t j = l; j < n; ++j) {
                        double s = 0.0;
                        for (size_t k = i; k < m; ++k) s += U[k][i] * U[k][j];
                        f = s / h;
                        for (size_t k = i; k < m; ++k) U[k][j] += f * U[k][i];
                    }
                    for (size_t k = i; k < m; ++k) U[k][i] *= scale;
                }
            }
            S[i] = scale * g;

            g = scale = 0.0;
            if (i < m && i != n - 1) {
                for (size_t k = l; k < n; ++k) scale += std::abs(U[i][k]);
                if (scale != 0.0) {
                    for (size_t k = l; k < n; ++k) {
                        U[i][k] /= scale;
                        g += U[i][k