机器学习算法SVD奇异值分解之c++实现（不调用外源库）

AI街潜水的八角

已于 2023-12-26 14:17:53 修改

阅读量450

点赞数

分类专栏：从零实现十大机器学习文章标签：机器学习算法 c++

于 2023-10-11 10:27:21 首次发布

本文链接：https://blog.csdn.net/m0_59023219/article/details/133762169

版权

从零实现十大机器学习专栏收录该内容

10 篇文章 2 订阅

订阅专栏

目前玩机器学习的小伙伴，上来就是使用现有的sklearn机器学习包，写两行代码，调调参数就能跑起来，看似方便，实则有时不利于个人能力发展，要知道现在公司需要的算法工程师，不仅仅只是会调参（这种工作，入门几个月的人就可以干了），而是要深入底层，能优化代码，能自己搭。

本文章适合以下几类人：

1）初学者，了解机器学习的实现过程

2）想提升自己的代码能力

第一步：原理

什么是SVD？

SVD(Singular Value Decomposition, 奇异值分解)是线性代数中既优雅又强大的工具, 它揭示了矩阵最本质的变换. 使用SVD对矩阵进行分解, 能得到代表矩阵最本质变化的矩阵元素. 这就好比一个合数能表示为若干质数之积, 分解合数能得到表示该合数的质因数; 复杂周期信号可以表示为若干简单的正弦波和余弦波之和, 使用傅里叶变换能得到表示该信号的简单波; 复杂矩阵所代表的线性变换可由若干个简单矩阵所代表的线性变换组合起来, 使用SVD能找到这些简单矩阵.

SVD的推导，百度一下很多，这里就做一下搬运工了，可参考：奇异值分解（SVD） - 知乎

第二步：代码实现

//还有一种奇异值分解，http://www.cnblogs.com/zhangchaoyang/articles/2575948.html

/**

函数原型:

bool svd(vector<vector<double> > A, int K, std::vector<std::vector<double> > &U, std::vector<double> &S, std::vector<std::vector<double> > &V);

输入矩阵A,分解矩阵的秩K

输出U,S,V

本函数将A分解为U diag(S) V'

S[i],U[i],V[i]是A的第i大奇异值，及其对应的左歧义向量和右奇异向量

S,U,V的size由K指定

K是需要分解的rank，0<K<=min(m,n)

本程序采用的是最基本幂迭代算法，在linux 下编译通过

**/

#include<math.h>
#include<iostream>
#include<string.h>
#include<vector>
#include<time.h>
#include <stdlib.h>
#include <string>
#include "matrix.h"
#include <fstream>
#include <sstream>
#include <stack>

using namespace std;
const int MAX_ITER=100000;
const double eps=0.0000001;

//累加，然后开方
double get_norm(double *x, int n)
{
    double r = 0;
    for (int i = 0; i < n; i++)
    {
        r += x[i] * x[i];
    }
    return sqrt(r);
}
//归一化
double normalize(double *x, int n)
{
    double r = get_norm(x,n);
    if (r < eps)
    {
        return 0;
    }
    for (int i = 0; i < n; i++)
    {
        x[i] /= r;
    }
    return r;
}
//两个不同的数相乘，然后累加
double product(double *a, double *b, int n)
{
    double r = 0;
    for (int i = 0; i < n; i++)
    {
        r += a[i] * b[i];
    }
    return r;
}

void orth(double *a, double *b, int n)
{
    double r = product(a,b,n);
    for (int i = 0; i < n; i++)
    {
        b[i] -= r*a[i];
    }
}
//A = U*S*VT U:m*m S:m*n V:n*n 正交矩阵 对角矩阵 正交矩阵
bool svd(vector<vector<double> > A, int K, vector<vector<double> > &U, vector<double> &S, vector<vector<double> > &V){

}

void print(vector<vector<double> > &A){

    for(int i=0;i<A.size();i++){

        for(int j=0;j<A[i].size();j++){

            cout<<A[i][j]<<' ';
        }
    cout<<endl;
    }
}

int main()
{
    dataToMatrix dtm;
    cout<<"loadData"<<endl;
    cout<<"----------------------"<<endl;
    char file[30]="G:/data/svd.txt";
    dtm.loadData(&dtm,file);
    dtm.print(dtm);
    Matrix x;
    x.loadMatrix(&x,dtm);

    int m=x.col;
    int n=x.row;

    srand(time(0));

    vector<vector<double> > A;
    A.resize(m);

    for(int i=0;i<m;i++){
        A[i].resize(n);
        for(int j=0;j<n;j++)
            A[i][j]=x.mat[i][j];
    }

    print(A);
    cout<<endl;

    vector<vector<double> > U;
    vector<double> S;
    vector<vector<double> > V;

    svd(A,16,U,S,V);
    cout<<"U="<<endl;

    print(U);

    cout<<endl;

    cout<<"S="<<endl;

    for(int i=0;i<S.size();i++){
    cout<<S[i]<<' ';
    }

    cout<<endl;

    cout<<"V="<<endl;

    print(V);

    //数组转换成矩阵
    Matrix UU;
    UU.initMatrix(&UU,32,2);
    Matrix VV;
    VV.initMatrix(&VV,2,32);
    Matrix SS;
    SS.initMatrix(&SS,2,2);
    Matrix result;
    result.initMatrix(&result,32,32);

    for (int i = 0; i < 32; i++)
    {
        for (int j = 0; j < 2; j++)
        {
            UU.mat[i][j] = U[j][i];
        }
    }

    SS.mat[0][0] = S[0];
    SS.mat[1][1] = S[1];

    for (int i = 0; i < 2; i++)
    {
        for (int j = 0; j < 32; j++)
        {
            VV.mat[i][j] = V[i][j];
        }
    }

    //将数据进行复原
    Matrix UUTemp;
    UUTemp.initMatrix(&UUTemp, 32, 2);

    UUTemp.multsmatrix(&UUTemp, UU,SS);
    result.multsmatrix(&result,UUTemp,VV);

    for (int i = 0; i < 32; i++)
    {
        for (int j = 0; j < 32; j++)
        {
            if (result.mat[i][j] > 0.5)//阈值
            {
                result.mat[i][j] = 1;
            }
            else
                result.mat[i][j] = 0;
        }
    }

    cout<<endl;
    result.print(result);

    return 0;

}

第三步：运行过程

运行结果