当前位置:   article > 正文

计算方法实验5:C++实现矩阵的奇异值分解

计算方法实验5:C++实现矩阵的奇异值分解

任务

生成一个4 × 3的随机矩阵A,应用Jacobi方法求解矩阵 A A T \mathbf{AA^T} AAT的特征值,计算矩阵A的SVD分解。要求A的每个元素均为[0:1]区间内的随机数。

算法

  1. Jacobi方法求矩阵 A A T \mathbf{AA^T} AAT的特征值;
  2. A A T \mathbf{AA^T} AAT的特征值从大到小排序后,用高斯消元法求解每个特征值对应的特征向量并归一化,然后转置为列向量组,得到 U 4 × 4 \mathbf{U}_{4\times4} U4×4;
  3. 矩阵 Σ 4 × 3 \mathbf{\Sigma}_{4\times3} Σ4×3的主对角元为 A A T \mathbf{AA^T} AAT的非零特征值的算术平方根(本题中仅考虑了实数情况),其余位置全为0;
  4. 矩阵 V T 3 × 3 \mathbf{V^T}_{3\times3} VT3×3每行即 U T \mathbf{U^T} UT A \mathbf{A} A对应行除以 Σ \mathbf{\Sigma} Σ对应的主对角元。

注:理论上 A A T \mathbf{AA^T} AAT A T A \mathbf{A^TA} ATA的非零特征值相同,但是程序计算出的结果会有一些差异,因此 Σ \mathbf{\Sigma} Σ的主对角元不能取 A T A \mathbf{A^TA} ATA的特征值的算术平方根, U T \mathbf{U^T} UT的每行元素也不能直接取 A T A \mathbf{A^TA} ATA的特征向量,而要按第4步所述计算,详情请参考这篇论文

代码

主函数

int main()
{
    vector<vector<long double>> A = generateRandomMatrix(4, 3);
    cout << "A: " << endl;
    for(int i = 0; i < 4; i++)
    {
        for(int j = 0; j < 3; j++)
            cout << A[i][j] << " ";
        cout << endl;
    }
    cout << endl;
    vector<vector<long double>> AT = calAT(A);
    vector<vector<long double>> AAT = multiplyMatrices(A, AT);
    int n1 = AAT.size();
    int n2 = A[0].size();
    vector<long double> x =calEigenValue(AAT);
    sort(x.begin(), x.end());
    reverse(x.begin(), x.end());

    vector<vector<long double>> U = calAT(Normalization(calEigenVector(AAT, x)));
    cout << "U: " << endl;
    for(int i = 0; i < n1; i++)
    {
        for(int j = 0; j < n1; j++)
            cout << U[i][j] << " ";
        cout << endl;
    }
    cout << endl;
    
    vector<vector<long double>> Sigma1 = calSigma(x, n1, n2);
    cout << "Sigma: " << endl;
    for(int i = 0; i < n1; i++)
    {
        for(int j = 0; j < n2; j++)
            cout << Sigma1[i][j] << " ";
        cout << endl;
    }
    cout << endl;

    vector<vector<long double>> VT = calculateVT(U, A, Sigma1, n2);
    cout << "VT: " << endl;
    for(int i = 0; i < n2; i++)
    {
        for(int j = 0; j < n2; j++)
            cout << VT[i][j] << " ";
        cout << endl;
    }
    return 0;
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49

生成一个具有指定行数和列数的随机矩阵

vector<vector<long double>> generateRandomMatrix(int rows, int cols) {
    random_device rd;//随机数种子
    mt19937 gen(rd());//使用 rd 生成的随机数来初始化一个mt19937类型的随机数生成器 gen
    uniform_real_distribution<> dis(0.0, 1.0);//创建一个均匀实数分布 dis

    vector<vector<long double>> matrix(rows, vector<long double>(cols));
    for (int i = 0; i < rows; ++i)
        for (int j = 0; j < cols; ++j)
            matrix[i][j] = dis(gen);//生成随机数
    return matrix;
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11

找到矩阵 A 中非对角元中绝对值最大的元素,并返回其位置

pair<int, int> chooseMax(vector<vector<long double>> A)
{
    long double max = 0;
    pair<int, int> maxPos;
    int n = A.size();
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++)
            if(i != j && fabsl(A[i][j]) > max)
            {
                max = fabsl(A[i][j]);
                maxPos = make_pair(i, j);
            }
    return maxPos;
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14

计算矩阵 A 的转置

vector<vector<long double>> calAT(vector<vector<long double>> A)
{
    int n1 = A.size();
    int n2 = A[0].size();
    vector<vector<long double>> AT(n2, vector<long double>(n1));
    for(int i = 0; i < n1; i++)
        for(int j = 0; j < n2; j++)
            AT[j][i] = A[i][j];
    return AT;
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10

计算两个矩阵的乘积

vector<vector<long double>> multiplyMatrices(const vector<vector<long double>> A, const vector<vector<long double>> B) {
    int n1 = A.size();
    int n2 = B[0].size();
    int n3 = A[0].size();
    vector<vector<long double>> result(n1, vector<long double>(n2, 0.0));

    for(int i = 0; i < n1; i++) {
        for(int j = 0; j < n2; j++) {
            for(int k = 0; k < n3; k++) {
                result[i][j] += A[i][k] * B[k][j];
            }
        }
    }
    return result;
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15

计算矩阵 Q^T * A * Q

vector<vector<long double>> calQTAQ(vector<vector<long double>> A)
{
    int n = A.size();
    pair<int, int> maxPos = chooseMax(A);
    int row = maxPos.first;
    int col = maxPos.second;
    
    long double s = (A[col][col] - A[row][row]) / (2 * A[row][col]);
    long double t = 0;
    if(s == 0)
        t = 1;
    else if(abs(-s + sqrt(1 + s * s)) <= abs(-s - sqrt(1 + s * s)))
        t = -s + sqrt(1 + s * s);
    else
        t = -s - sqrt(1 + s * s);

    long double c = 1 / sqrt(1 + t * t);
    long double d = t * c;

    vector<vector<long double>> QTAQ(n, vector<long double>(n));
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++)
            QTAQ[i][j] = calculateElement(A, i, j, row, col, t, c, d);
    return QTAQ;
}

// 计算矩阵Q^T * A * Q的每个元素,使用给定的参数 p, q, t, c, d
long double calculateElement(const vector<vector<long double>> A, int i, int j, long double p, long double q, long double t, long double c, long double d) {
    if (i == p && j == p)
        return A[p][p] - t * A[p][q];
    else if (i == q && j == q)
        return A[q][q] + t * A[p][q];
    else if ((i == p && j == q) || (i == q && j == p))
        return 0;
    else if (i != q && i != p && (j == p || j == q))
        return (j == p ? c : d) * A[p][i] - (j == p ? d : (-c)) * A[q][i];
    else if ((i == p || i == q) && j != q && j != p)
        return (i == p ? c : d) * A[p][j] - (i == p ? d : (-c)) * A[q][j];
    else
        return A[i][j];
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41

判断Jacobi方法是否满足结束条件

int judgeEnd(vector<vector<long double>> A)
{
    int i, j;
    int n = A.size();
    for(i = 0; i < n; i++)
        for(j = 0; j < n; j++)
            if(i != j && fabsl(A[i][j]) >= 1e-6)
                return 0;
    if(i == n && j == n) 
        return 1;
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11

Jacobi方法计算矩阵 A 的特征值

vector<long double> calEigenValue(vector<vector<long double>> A)
{
    int n = A.size();
    vector<long double> eigenValue(n);
    vector<vector<long double>> QTAQ= calQTAQ(A);
    int i, j;
    while(!judgeEnd(QTAQ))
        QTAQ = calQTAQ(QTAQ);
    for(i = 0; i < n; i++)
        eigenValue[i] =QTAQ[i][i];
    return eigenValue;
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12

使用给定的特征值计算矩阵 A 的特征向量

vector<vector<long double>> calEigenVector(vector<vector<long double>> A, vector<long double> eigenValue)
{
    int n = A.size();
    int num = 0;
    vector<vector<long double>> x(n, vector<long double>(n));
    vector<vector<long double>> tempMartix(n, vector<long double>(n));
    vector<vector<long double>> eigenVector(n, vector<long double>(n));
    for(int k = 0; k < n; k++)
    {
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++)
                i == j ? tempMartix[i][j] = A[i][j] - eigenValue[k] : tempMartix[i][j] = A[i][j];

        vector<vector<long double>> B = Column_Elimination(tempMartix);
        int cnt = 0;//记录消元后全为0的行数
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
            {
                if(fabsl(B[i][j]) > 1e-7)
                    break;
                else if(j == n - 1)
                    cnt++;
            }
        }
        vector<vector<long double>> result = solve(B, cnt);
        for(int i = 0; i < cnt; i++)
            copy(result[i].begin(), result[i].end(), x[num + i].begin());
        num += cnt;
    }
    return x;
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32

求解有无穷多解的线性方程组

// 对矩阵A进行列主元化成上三角
vector<vector<long double>> Column_Elimination(vector<vector<long double>> A)
{
    int n = A.size();
    vector<vector<long double>> Temp(n, vector<long double>(n));
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++)
            Temp[i][j] = A[i][j];
    for(int col = 0; col < n; col++)
    {
        long double maxnum = abs(Temp[col][col]);
        int maxrow = col;
        for(int row = col + 1; row < n; row++)
        {
            if(abs(Temp[row][col]) > maxnum)
            {
                maxnum = abs(Temp[row][col]);
                maxrow = row;
            }
        }
        swap(Temp[col], Temp[maxrow]);
        for(int row = col + 1; row < n; row++)
        {
            long double res = Temp[row][col] / Temp[col][col];
            for(int loc = col; loc < n; loc++)
                Temp[row][loc] -= Temp[col][loc] * res; 
        }
    }
    return Temp;
}

// 求解系数矩阵为上三角矩阵A,且A的行列式不为0的线性方程组
vector<long double> SolveUpperTriangle(vector<vector<long double>> A, vector<long double> b)
{
    int n = A.size();
    vector<long double> x(n);
    vector<vector<long double>> Temp(n, vector<long double>(n+1));
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++)
            Temp[i][j] = A[i][j];
    for(int i = 0; i < n; i++)
        Temp[i][n] = b[i];
    for(int row = n-1; row >= 0; row--)
    {
        for(int col = row + 1; col < n; col++)
        {
            Temp[row][n] -= Temp[col][n] * Temp[row][col] / Temp[col][col];
            Temp[row][col] = 0;
        }
        Temp[row][n] /= Temp[row][row];
        Temp[row][row] = 1;
    }
    for(int i = 0; i < n; i++)
        x[i] = Temp[i][n];
    return x;
}

// 解系数矩阵为上三角矩阵 A 的线性方程组,且A全为0的行数为 cnt
vector<vector<long double>> solve(vector<vector<long double>> A, int cnt)
{
    int n = A.size();
    vector<vector<long double>> x(cnt, vector<long double>(n));
    vector<vector<long double>> Temp(n-cnt, vector<long double>(n-cnt));
    vector<long double> Tempb(n-cnt);
    for(int i = 0; i < cnt; i++)
    {
        for(int j = n - 1; j >= n - cnt; j--)
        {
            if(j >= n - i)
                x[i][j] = 0;
            else
                x[i][j] = 1;
        }
    }
    for(int i = 0; i < n - cnt; i++)
        for(int j = 0; j < n - cnt; j++)
            Temp[i][j] = A[i][j];
    for(int i = 0; i < cnt; i++)
    {
        for(int j = n - cnt - 1; j >=  0; j--)
        {
            Tempb[j] = 0;
            for(int k = 0; k < cnt; k++)
                Tempb[j] -= A[j][n- cnt + k] * x[i][n- cnt + k];
        }
        vector<long double> res = SolveUpperTriangle(Temp, Tempb);
        for(int j = 0; j < n - cnt; j++)
            x[i][j] = res[j];
    }
    return x;
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76
  • 77
  • 78
  • 79
  • 80
  • 81
  • 82
  • 83
  • 84
  • 85
  • 86
  • 87
  • 88
  • 89
  • 90
  • 91

使用给定的特征值 x 和矩阵的行数 n1 和列数 n2,计算 Sigma 矩阵

vector<vector<long double>> calSigma(vector<long double> x, int n1, int n2)
{
    vector<vector<long double>> Sigma(n1, vector<long double>(n2));
    for(int i = 0; i < min(n1, n2); i++)
        Sigma[i][i] = sqrt(x[i]);
    return Sigma;
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7

计算向量 x 的欧几里得范数

long double EuclideanNorm(vector<long double> x)
{
    long double norm = 0;
    for(int i = 0; i < x.size(); i++)
        norm += x[i] * x[i];
    return sqrt(norm);
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7

对矩阵 A 进行归一化

vector<vector<long double>> Normalization(vector<vector<long double>> A)
{
    int rows = A.size();
    for(int i = 0; i < rows; i++)
    {
        long double norm = EuclideanNorm(A[i]);
        int cols = A[i].size();
        for(int j = 0; j < cols; j++)
            A[i][j] /= norm;
    }
    return A;
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12

计算 VT 矩阵

vector<vector<long double>> calculateVT(const vector<vector<long double>>& U, const vector<vector<long double>>& A, const vector<vector<long double>>& Sigma1, int n2) {
    vector<vector<long double>> UTA = multiplyMatrices(calAT(U), A);
    vector<vector<long double>> VT(n2, vector<long double>(n2));
    for(int i = 0; i < n2; i++)
        for(int j = 0; j < n2; j++)
            VT[i][j] = UTA[i][j]/Sigma1[i][i];
    return VT;
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8

附录

1.python实现特征值分解代码

import numpy as np
import math 

# Define the A2 matrix
A2 = np.array([
    [1.07251, 1.58415, 0.811742],
    [1.58415, 2.41048, 1.24756],
    [0.811742, 1.24756, 0.796375],
])

# Compute the eigenvalues of A2
eigenvalues_A2, eigenvectors_A2 = np.linalg.eig(A2)

print("sqrt of Eigenvalues of A2:", np.sqrt(eigenvalues_A2))
print("Eigenvectors of A2:")
for i in range(len(eigenvalues_A2)):
    print("Eigenvector for eigenvalue {}: {}".format(eigenvalues_A2[i], eigenvectors_A2[:, i]))
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17

2.python实现SVD分解代码

import numpy as np
from scipy.linalg import svd

# 定义矩阵
A = np.array([[0.726249, 0.915339, 0.594324],
              [0.515358, 0.975149, 0.661561],
              [0.528652, 0.788493, 0.0741007],
              [0.32985, 0.583744, 0.0309722]])

# 进行奇异值分解
U, s, VT = svd(A)

print("U:\n", U)
print("s:\n", s)
print("VT:\n", VT)
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15

3.C++判断一个向量是否为某个特征值对应的特征向量

#include<bits/stdc++.h>
using namespace std;
bool isEigenVector(vector<vector<long double>> A, vector<long double> v, long double lambda) {
    int n = A.size();
    vector<long double> Av(n);
    for(int i = 0; i < n; i++) {
        for(int j = 0; j < n; j++) {
            Av[i] += A[i][j] * v[j];
        }
        if(abs(Av[i] - lambda * v[i]) > 1e-4) {
            return false;
        }
    }
    return true;
}

int main()
{
    vector<vector<long double>> A = {{0.68112222, -0.03900667, 1.26519111, 0.51345778},
                                     {-0.03900667, 0.18675067, -0.319568, -0.11719467},
                                     {1.26519111, -0.319568, 3.09242489, 1.28774489},
                                     {0.51345778, -0.11719467, 1.28774489, 0.57853156}};
    
    vector<long double> v = {-0.47971899,  0.07252408 , 0.1757674,0.85657211};
    long double lambda = 4.196675163197978;
    if(isEigenVector(A, v, lambda))
        cout << "v is an eigenvector of A with eigenvalue lambda." << endl;
    else
        cout << "v is not an eigenvector of A with eigenvalue lambda." << endl;
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/菜鸟追梦旅行/article/detail/458978
推荐阅读
相关标签
  

闽ICP备14008679号