计算机视觉
图像处理

机器学习—奇异值分解

讲解主要内容如下:

1.什么是奇异值分解?为什么任意实矩阵都存在奇异值分解?

2.怎么用C语言代码实现SVD分解?

3.实际应用:

基于SVD的图像压缩

基于SVD的协同过滤推荐系统

一、SVD奇异值分解概念

在多数情况下,数据中的一小段携带了数据集中的大部分信息,其他信息要么是噪声,要么就是毫不相干的信息。在线性代数中有很多矩阵分解技术,通过它们可以将原始矩阵表示成新的易于处理的形式。不同的矩阵分解技术具有不同的性质,有各自特定的应用场景。

奇异值分解SVD作为矩阵分解的一种类型,可以使我们只用很小的数据集就能表达原始数据集,它可以从噪声数据中抽取特征值。SVD可以将任意维数的原始数据集矩阵Data分解成三个矩阵U、Σ、VT,如下公式所示:

上述分解中会构建出一个矩阵Σ,该矩阵只有对角元素,其他元素均为0.此外,Σ的对角元素是从大到小排列的,这些对角线元素称为原始数据集矩阵Data的 “奇异值”singular value。它与PCA特征的联系:PCA得到的是矩阵的特征值,奇异值就是矩阵Data*DataT特征值的平方根。

设A∈Rrm*n,ATA的特征值为:


则称  (i=1,2,…,n)为矩阵A的奇异值。

当A为零矩阵时,它的所有奇异值均为零。一般的,矩阵A的奇异值个数等于A的列数,A的非零奇异值的个数等于A的秩。

接下来,首先用数学推导介绍SVD存在的必然性。

一些定义:

正交矩阵

如果方正Q满足QTQ=QQT=I(或者等价的Q-1=QT),则称Q为正交矩阵。

非奇异矩阵

若n阶矩阵A的行列式不为零,即 |A|≠0,则称A为非奇异矩阵或满秩矩阵,否则称A为奇异矩阵或降秩矩阵。

正定矩阵

设M是n阶方阵,如果对任何非零向量z,都有 zTMz > 0,就称M正定矩阵。

定理证明:

定理1 (实对称矩阵正交相似于对角矩阵)

对任意n阶实对称矩阵A,都有n阶正交矩阵Q,使得


为了证明该定理,首先给出并证明两个引理:

  • . 实对称矩阵特征值为实数
  • 若任意n阶矩阵的特征值为实数,则有正交矩阵Q,使得

引理1,证明:

设λ为A的一个特征值,x为对应的特征向量,则有Ax = λx,因此

所以特征值必为实数。

引理2,利用数学归纳法证明:

当n=1时,结论显然成立

设对n-1阶矩阵的结论也成立,现在证明对n阶矩阵结论依然成立

令λ1为A的一个实特征值,相应的特征向量为a1,不妨设a1已单位化。把a1扩充为Rn的标准正交基a1, a2,…, an,构造矩阵X=(a2,…, an)和P=(a1,X),则P为正交矩阵且有:

于是,AP = (Aa1, AX)= (λ1a1, P(P-1AX)) = (λ11,P(P-1AX)) = P(λ1ε1, P-1AX)。设:

个人理解:A->n*n,X->n*n-1,P->n*n,则P-1AX->n*n-1,正好这么拆,a->1*n-1

从而有:

根据归纳假设,对于A1有n-1阶正交矩阵T,使得T-1A1T = TTA1T= At为上三角矩阵,现取:

由于P和Q都是正交矩阵,而正交矩阵相乘结果R仍是正交矩阵,因此可以写作:


现在回到定理1的证明,因为A为实对称矩阵,则有引理1知,A的特征值均为实数。由引理2知,存在正交矩阵R,使得

C为上三角矩阵,而A=AT,则有C=CT,但C是上三角矩阵,CT为下三角矩阵,所以C必为对角矩阵。定理1得证。其中,C的对角线元素即为A的特征值构成。

定理2(奇异值分解)

设A∈Rrm*n,则存在m阶正交矩阵U和n阶正交矩阵V,使得

其中,Σ = diag(σ12,…, σr),即σ1, σ2,…,σr是A的非零奇异值,亦记作A=UDVT

证明:

记对称矩阵ATA的特征值为:

由对称矩阵的正交分解(定理1),我们有:


由第一个公式可以得到,

由第二个公式可以得到,


故有,证毕

若记U=(u1, u2,…,um),V=(v1, v2,…,vn),根据SVD的分解,A可以表示为:

上式就是A的奇异值分解过程。

【我恨死CSDN的公式编辑了!!每次都贴图,烦的要命】

二、SVD奇异值分解求解及实现

上面的一堆推导,只是为了说明对任意矩阵都存在奇异值分解,而没有讲述如何求解奇异值分解。以下内容来自《徐士良C常用算法程序集(第二版)》,一般实矩阵的奇异值分解:

用household变换和变形QR算法对一般实矩阵进行奇异值分解,得到左右奇异向量U、V及奇异值矩阵D。

[个人认为具体怎么求SVD,大家不要去记和理解了,要用的时候能找到源码,方便其他移植即可]

例:求下面两个矩阵A和B的奇异值分解,ε取0.000001

C语言代码:

徐士良老头写的,可读性很差

bmuav.c:

#include "stdlib.h"
#include "math.h"
void ppp(double *a, double *e,double *s,double *v,int m,int n);
void sss(double *fg,double *cs);
/************************************************************************/
/* input:
/* a:存放m*n实矩阵A,返回时亦是奇异矩阵
/* m:行数 n:列数
/* u:存放m*m左奇异向量, v:存放n*n右奇异向量 
/* eps:给定精度要求,  ka: max(m,n)+1
/* output:
/* 返回值如果为负数,表示迭代了60次,还未求出奇异值;返回值为非负数,正常
/************************************************************************/
int bmuav(double *a,int m,int n,double *u,double *v,double eps,int ka)
{ 
	int i,j,k,l,it,ll,kk,ix,iy,mm,nn,iz,m1,ks;
	double d,dd,t,sm,sm1,em1,sk,ek,b,c,shh,fg[2],cs[2];
	double *s,*e,*w;

	s = (double *)malloc(ka*sizeof(double));
	e = (double *)malloc(ka*sizeof(double));
	w = (double *)malloc(ka*sizeof(double));
	it=60; 
	k=n;
	
	if (m-1<n) 
		k=m-1;
	
	l=m;

	if (n-2<m)
		l=n-2;

	if (l<0) l=0; ll=k; if (l>k) 
		ll=l;

	if (ll>=1)
	{ 
		for (kk=1; kk<=ll; kk++)
		{	
			if (kk<=k)
			{	d=0.0;
				for (i=kk; i<=m; i++)
				{
					ix=(i-1)*n+kk-1; d=d+a[ix]*a[ix];
				}
				s[kk-1]=sqrt(d);
				if (s[kk-1]!=0.0)
				{
					ix=(kk-1)*n+kk-1;
					if (a[ix]!=0.0)
					{
						s[kk-1]=fabs(s[kk-1]);
						if (a[ix]<0.0) 
							s[kk-1]=-s[kk-1];
					}
					for (i=kk; i<=m; i++) { iy=(i-1)*n+kk-1; a[iy]=a[iy]/s[kk-1]; } a[ix]=1.0+a[ix]; } s[kk-1]=-s[kk-1]; } if (n>=kk+1)
			{ 
				for (j=kk+1; j<=n; j++)
				{
					if ((kk<=k)&&(s[kk-1]!=0.0))
					{
						d=0.0;
						for (i=kk; i<=m; i++)
						{ 
							ix=(i-1)*n+kk-1;
							iy=(i-1)*n+j-1;
							d=d+a[ix]*a[iy];
						}
						d=-d/a[(kk-1)*n+kk-1];
						for (i=kk; i<=m; i++)
						{ 
							ix=(i-1)*n+j-1;
							iy=(i-1)*n+kk-1;
							a[ix]=a[ix]+d*a[iy];
						}
					}
					e[j-1]=a[(kk-1)*n+j-1];
				}
			}
			if (kk<=k)
			{ 
				for (i=kk; i<=m; i++)
				{ 
					ix=(i-1)*m+kk-1; iy=(i-1)*n+kk-1;
					u[ix]=a[iy];
				}
			}
			if (kk<=l)
			{ 
				d=0.0;
				for (i=kk+1; i<=n; i++)
					d=d+e[i-1]*e[i-1];
				e[kk-1]=sqrt(d);
				if (e[kk-1]!=0.0)
				{
					if (e[kk]!=0.0)
					{ 
						e[kk-1]=fabs(e[kk-1]);
						if (e[kk]<0.0)
							e[kk-1]=-e[kk-1];
					}
					for (i=kk+1; i<=n; i++)
						e[i-1]=e[i-1]/e[kk-1];
					e[kk]=1.0+e[kk];
				}
				e[kk-1]=-e[kk-1];
				if ((kk+1<=m)&&(e[kk-1]!=0.0))
				{ 
					for (i=kk+1; i<=m; i++) w[i-1]=0.0;
						for (j=kk+1; j<=n; j++)
							for (i=kk+1; i<=m; i++)
								w[i-1]=w[i-1]+e[j-1]*a[(i-1)*n+j-1];
					for (j=kk+1; j<=n; j++)
						for (i=kk+1; i<=m; i++)
						{ 
							ix=(i-1)*n+j-1;
							a[ix]=a[ix]-w[i-1]*e[j-1]/e[kk];
						}
				}
				for (i=kk+1; i<=n; i++)
					v[(i-1)*n+kk-1]=e[i-1];
			}
		}
	}
	mm=n;
	if (m+1<n) 
		mm=m+1;
	if (k<n) 
		s[k]=a[k*n+k];
	if (m<mm) 
		s[mm-1]=0.0;
	if (l+1<mm) e[l]=a[l*n+mm-1]; e[mm-1]=0.0; nn=m; if (m>n) 
		nn=n;
	if (nn>=k+1)
	{ 
		for (j=k+1; j<=nn; j++)
		{ 
			for (i=1; i<=m; i++) u[(i-1)*m+j-1]=0.0; u[(j-1)*m+j-1]=1.0; } } if (k>=1)
	{ 
		for (ll=1; ll<=k; ll++) { kk=k-ll+1; iz=(kk-1)*m+kk-1; if (s[kk-1]!=0.0) { if (nn>=kk+1)
				for (j=kk+1; j<=nn; j++)
				{
					d=0.0;
					for (i=kk; i<=m; i++)
					{ 
						ix=(i-1)*m+kk-1;
						iy=(i-1)*m+j-1;
						d=d+u[ix]*u[iy]/u[iz];
					}
					d=-d;
					for (i=kk; i<=m; i++)
					{ 
						ix=(i-1)*m+j-1;
						iy=(i-1)*m+kk-1;
						u[ix]=u[ix]+d*u[iy];
					}
				}
				for (i=kk; i<=m; i++) { ix=(i-1)*m+kk-1; u[ix]=-u[ix]; } u[iz]=1.0+u[iz]; if (kk-1>=1)
					for (i=1; i<=kk-1; i++)
						u[(i-1)*m+kk-1]=0.0;
			}
			else
			{ 
				for (i=1; i<=m; i++)
					u[(i-1)*m+kk-1]=0.0;
				u[(kk-1)*m+kk-1]=1.0;
			}
		}
	}
	for (ll=1; ll<=n; ll++)
	{ 
		kk=n-ll+1; iz=kk*n+kk-1;
		if ((kk<=l)&&(e[kk-1]!=0.0))
		{ 
			for (j=kk+1; j<=n; j++)
			{
				d=0.0;
				for (i=kk+1; i<=n; i++)
				{ 
					ix=(i-1)*n+kk-1; iy=(i-1)*n+j-1;
					d=d+v[ix]*v[iy]/v[iz];
				}
				d=-d;
				for (i=kk+1; i<=n; i++)
				{ 
					ix=(i-1)*n+j-1; 
					iy=(i-1)*n+kk-1;
					v[ix]=v[ix]+d*v[iy];
				}
			}
		}
		for (i=1; i<=n; i++)
			v[(i-1)*n+kk-1]=0.0;
		v[iz-n]=1.0;
	}
	for (i=1; i<=m; i++)
		for (j=1; j<=n; j++) a[(i-1)*n+j-1]=0.0; m1=mm; it=60; while (1==1) { if (mm==0) { ppp(a,e,s,v,m,n); free(s); free(e); free(w); return(1); } if (it==0) { ppp(a,e,s,v,m,n); free(s); free(e); free(w); return(-1); } kk=mm-1; while ((kk!=0)&&(fabs(e[kk-1])!=0.0)) { d=fabs(s[kk-1])+fabs(s[kk]); dd=fabs(e[kk-1]); if (dd>eps*d) 
				kk=kk-1;
			else
				e[kk-1]=0.0;
		}
		if (kk==mm-1)
		{ 
			kk=kk+1;
			if (s[kk-1]<0.0)
			{
				s[kk-1]=-s[kk-1];
				for (i=1; i<=n; i++)
				{
					ix=(i-1)*n+kk-1; v[ix]=-v[ix];
				}
			}
			while ((kk!=m1)&&(s[kk-1]<s[kk]))
			{ 
				d=s[kk-1]; s[kk-1]=s[kk]; s[kk]=d;
				if (kk<n)
					for (i=1; i<=n; i++)
					{
						ix=(i-1)*n+kk-1; iy=(i-1)*n+kk;
						d=v[ix]; v[ix]=v[iy]; v[iy]=d;
					}
				if (kk<m)
				for (i=1; i<=m; i++) { ix=(i-1)*m+kk-1; iy=(i-1)*m+kk; d=u[ix]; u[ix]=u[iy]; u[iy]=d; } kk=kk+1; } it=60; mm=mm-1; } else { ks=mm; while ((ks>kk)&&(fabs(s[ks-1])!=0.0))
			{
				d=0.0;
				if (ks!=mm)
					d=d+fabs(e[ks-1]);
				if (ks!=kk+1)
					d=d+fabs(e[ks-2]);
				dd=fabs(s[ks-1]);
				if (dd>eps*d)
					ks=ks-1;
				else
					s[ks-1]=0.0;
			}
			if (ks==kk)
			{
				kk=kk+1;
				d=fabs(s[mm-1]);
				t=fabs(s[mm-2]);
				if (t>d)
					d=t;
				t=fabs(e[mm-2]);
				if (t>d)
					d=t;
				t=fabs(s[kk-1]);
				if (t>d)
					d=t;
				t=fabs(e[kk-1]);
				if (t>d)
					d=t;
				sm=s[mm-1]/d;
				sm1=s[mm-2]/d;
				em1=e[mm-2]/d;
				sk=s[kk-1]/d;
				ek=e[kk-1]/d;
				b=((sm1+sm)*(sm1-sm)+em1*em1)/2.0;
				c=sm*em1; c=c*c; shh=0.0;
				if ((b!=0.0)||(c!=0.0))
				{ 
					shh=sqrt(b*b+c);
					if (b<0.0)
						shh=-shh;
					shh=c/(b+shh);
				}
				fg[0]=(sk+sm)*(sk-sm)-shh;
				fg[1]=sk*ek;
				for (i=kk; i<=mm-1; i++)
				{ 
					sss(fg,cs);
					if (i!=kk)
						e[i-2]=fg[0];
					fg[0]=cs[0]*s[i-1]+cs[1]*e[i-1];
					e[i-1]=cs[0]*e[i-1]-cs[1]*s[i-1];
					fg[1]=cs[1]*s[i];
					s[i]=cs[0]*s[i];
					if ((cs[0]!=1.0)||(cs[1]!=0.0))
						for (j=1; j<=n; j++)
						{
							ix=(j-1)*n+i-1;
							iy=(j-1)*n+i;
							d=cs[0]*v[ix]+cs[1]*v[iy];
							v[iy]=-cs[1]*v[ix]+cs[0]*v[iy];
							v[ix]=d;
						}
					sss(fg,cs);
					s[i-1]=fg[0];
					fg[0]=cs[0]*e[i-1]+cs[1]*s[i];
					s[i]=-cs[1]*e[i-1]+cs[0]*s[i];
					fg[1]=cs[1]*e[i];
					e[i]=cs[0]*e[i];
					if (i<m)
						if ((cs[0]!=1.0)||(cs[1]!=0.0))
							for (j=1; j<=m; j++)
							{ 
								ix=(j-1)*m+i-1;
								iy=(j-1)*m+i;
								d=cs[0]*u[ix]+cs[1]*u[iy];
								u[iy]=-cs[1]*u[ix]+cs[0]*u[iy];
								u[ix]=d;
							}
				}
				e[mm-2]=fg[0];
				it=it-1;
			}
			else
			{ 
				if (ks==mm)
				{ 
					kk=kk+1;
					fg[1]=e[mm-2]; e[mm-2]=0.0;
					for (ll=kk; ll<=mm-1; ll++)
					{ 
						i=mm+kk-ll-1;
						fg[0]=s[i-1];
						sss(fg,cs);
						s[i-1]=fg[0];
						if (i!=kk)
						{
							fg[1]=-cs[1]*e[i-2];
							e[i-2]=cs[0]*e[i-2];
						}
						if ((cs[0]!=1.0)||(cs[1]!=0.0))
							for (j=1; j<=n; j++)
							{ 
								ix=(j-1)*n+i-1;
								iy=(j-1)*n+mm-1;
								d=cs[0]*v[ix]+cs[1]*v[iy];
								v[iy]=-cs[1]*v[ix]+cs[0]*v[iy];
								v[ix]=d;
							}
					}
				}
				else
				{ 
					kk=ks+1;
					fg[1]=e[kk-2];
					e[kk-2]=0.0;
					for (i=kk; i<=mm; i++)
					{ 
						fg[0]=s[i-1];
						sss(fg,cs);
						s[i-1]=fg[0];
						fg[1]=-cs[1]*e[i-1];
						e[i-1]=cs[0]*e[i-1];
						if ((cs[0]!=1.0)||(cs[1]!=0.0))
							for (j=1; j<=m; j++) { ix=(j-1)*m+i-1; iy=(j-1)*m+kk-2; d=cs[0]*u[ix]+cs[1]*u[iy]; u[iy]=-cs[1]*u[ix]+cs[0]*u[iy]; u[ix]=d; } } } } } } return(1); } static void ppp(double a[], double e[],double s[],double v[],int m,int n) { int i,j,p,q; double d; if (m>=n) i=n;
	else i=m;
	for (j=1; j<=i-1; j++)
	{
		a[(j-1)*n+j-1]=s[j-1];
		a[(j-1)*n+j]=e[j-1];
	}
	a[(i-1)*n+i-1]=s[i-1];
	if (m<n) 
		a[(i-1)*n+i]=e[i-1];
	for (i=1; i<=n-1; i++)
		for (j=i+1; j<=n; j++) { p=(i-1)*n+j-1; q=(j-1)*n+i-1; d=v[p]; v[p]=v[q]; v[q]=d; } return; } static void sss(double fg[],double cs[]) { double r,d; if ((fabs(fg[0])+fabs(fg[1]))==0.0) { cs[0]=1.0; cs[1]=0.0; d=0.0; } else { d=sqrt(fg[0]*fg[0]+fg[1]*fg[1]); if (fabs(fg[0])>fabs(fg[1]))
		{ 
			d=fabs(d);
			if (fg[0]<0.0) d=-d; } if (fabs(fg[1])>=fabs(fg[0]))
		{ 
			d=fabs(d);
			if (fg[1]<0.0) d=-d; } cs[0]=fg[0]/d; cs[1]=fg[1]/d; } r=1.0; if (fabs(fg[0])>fabs(fg[1])) 
		r=cs[1];
	else
		if (cs[0]!=0.0) r=1.0/cs[0];
	fg[0]=d;
	fg[1]=r;
	return;
}

brmul.c:矩阵乘法

调用:

// svd.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include "bmuav.c"
#include "brmul.c"

extern void brmul(double *a,double *b,int m,int n,int k,double *c);
extern int bmuav(double *a,int m,int n,double *u,double *v,double eps,int ka);

int _tmain(int argc, _TCHAR* argv[])
{
	int i,j;
    static double a[4][3]={ {1.0,1.0,-1.0},{2.0,1.0,0.0},
                           {1.0,-1.0,0.0},{-1.0,2.0,1.0}};
    static double b[3][4]={ {1.0,1.0,-1.0,-1.0},{2.0,1.0,
                            0.0,2.0},{1.0,-1.0,0.0,1.0}};
    static double u[4][4],v[3][3],c[4][3],d[3][4];
    double eps;
    eps=0.000001;
    i=bmuav(&a[0][0],4,3,&u[0][0],&v[0][0],eps,5);
    printf("\n");
    printf("EXAMPLE(1)\n");
    printf("\n");
    printf("i=%d\n",i);
    printf("\n");
    printf("MAT U IS:\n");
    for (i=0; i<=3; i++)
      { for (j=0; j<=3; j++)
          printf("%13.7e ",u[i][j]);
        printf("\n");
      }
    printf("\n");
    printf("MAT V IS:\n");
    for (i=0; i<=2; i++)
      { for (j=0; j<=2; j++)
          printf("%13.7e ",v[i][j]);
        printf("\n");
      }
    printf("\n");
    printf("MAT A IS:\n");
    for (i=0; i<=3; i++)
      { for (j=0; j<=2; j++)
          printf("%13.7e ",a[i][j]);
        printf("\n");
      }
    printf("\n\n");
    printf("MAT UAV IS:\n");
    brmul(&u[0][0],&a[0][0],4,4,3,&c[0][0]);
    brmul(&c[0][0],&v[0][0],4,3,3,&a[0][0]);
    for (i=0; i<=3; i++)
      { for (j=0; j<=2; j++)
          printf("%13.7e ",a[i][j]);
        printf("\n");
      }
    printf("\n\n");
    printf("EXAMPLE(2)\n");
    printf("\n");
    i=bmuav(&b[0][0],3,4,&v[0][0],&u[0][0],eps,5);
    printf("i=%d\n",i);
    printf("\n");
    printf("MAT U IS:\n");
    for (i=0; i<=2; i++)
      { for (j=0; j<=2; j++)
          printf("%13.7e ",v[i][j]);
        printf("\n");
      }
    printf("\n");
    printf("MAT V IS:\n");
    for (i=0; i<=3; i++)
      { for (j=0; j<=3; j++)
          printf("%13.7e ",u[i][j]);
        printf("\n");
      }
    printf("\n");
    printf("MAT B IS:\n");
    for (i=0; i<=2; i++)
      { for (j=0; j<=3; j++)
          printf("%13.7e ",b[i][j]);
        printf("\n");
      }
    printf("\n\n");
    printf("MAT UBV IS:\n");
    brmul(&v[0][0],&b[0][0],3,3,4,&d[0][0]);
    brmul(&d[0][0],&u[0][0],3,4,4,&b[0][0]);
    for (i=0; i<=2; i++)
      { for (j=0; j<=3; j++)
          printf("%13.7e ",b[i][j]);
        printf("\n");
      }
    printf("\n");
	return 0;
}

程序结果:

三、SVD应用实例

1.      基于SVD的图像压缩

这个例子比较简单,首先进行奇异值分解,得到奇异值矩阵,和左右奇异向量。然后由于只要很少的奇异值,就能包含绝大部分被分解的矩阵信息,因此我们挑选不 同数量的奇异值,重构图像,比较差异。这边分别实现了灰度图、RGB三色图的SVD分解。【奇异值到底选多少,自己打印奇异值矩阵,从大到小排序的,小到什么程度就舍弃,实际情况实际操作。。。】

#include "stdafx.h"
#include "cv.h"
#include "highgui.h"
#include "bmuav.c"
#include "brmul.c"
#define max(a,b)            (((a) > (b)) ? (a) : (b))

extern void brmul(double *a,double *b,int m,int n,int k,double *c);
extern int bmuav(double *a,int m,int n,double *u,double *v,double eps,int ka);
int Process(IplImage *src);

int _tmain(int argc, _TCHAR* argv[])
{
	IplImage *src = cvLoadImage("test1.jpg", CV_LOAD_IMAGE_GRAYSCALE);
	Process(src);
	return 0;
}

int Process(IplImage *src)
{
	double *data, *u, *v, *c;
	int height, width, i, j, ka;
	double eps;
	int scale;
	IplImage *dst;
	
	dst = cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,1);

	eps=0.000001;
	height = src->height;
	width = src->widthStep;
	
	//allocate memory for matrix
	data = (double *)malloc(sizeof(double)*height*width);
	u = (double *)malloc(sizeof(double)*height*height);
	v = (double *)malloc(sizeof(double)*width*width);
	memset(u,0,sizeof(double)*height*height);
	memset(v,0,sizeof(double)*width*width);
	
	if(NULL == data || NULL == u || NULL == v)
	{
		return -1;
	}

	//assign value
	for(i = 0;i < height;i++)
	{
		for(j = 0;j < width;j++) { data[i*width+j] = (double)(unsigned char)src->imageData[i*width+j];
		}
	}
	
	ka = max(height,width) + 1;
	bmuav(data, height, width, u, v, eps, ka);
	
	//dump svd, scale is selected by watching top xxx large data
	/*for (i=0; i<=100; i++)
	{ 
		for (j=0; j<=100; j++)
			printf("%f", data[i*width+j]);
		printf("\n");
	}*/
	
	//reconstruction
	scale = 50;
	for(i = scale;i<height;i++)
	{
		data[i*width+i] = 0;
	}

	/*c needs to be initilized here ,but in matrix mutiply funciton*/
	c = (double *)malloc(sizeof(double)*height*width);
	
	brmul(u, data ,height, height, width, c);
	brmul(c, v, height, width, width, data);

	//assign value
	for(i = 0;i < height;i++)
	{
		for(j = 0;j < width;j++) { dst->imageData[i*width+j] = (unsigned char)data[i*width+j];
		}
	}
	cvSaveImage("result.jpg",dst);
	free(data);free(u);free(v);
	cvReleaseImage(&dst);
	return 0;
}

结果:

彩色图:

#include "stdafx.h"
#include "cv.h"
#include "highgui.h"
#include "bmuav.c"
#include "brmul.c"
#define max(a,b)            (((a) > (b)) ? (a) : (b))

extern void brmul(double *a,double *b,int m,int n,int k,double *c);
extern int bmuav(double *a,int m,int n,double *u,double *v,double eps,int ka);
int Process(IplImage *src);

int _tmain(int argc, _TCHAR* argv[])
{
	IplImage *src = cvLoadImage("test3.jpg", CV_LOAD_IMAGE_UNCHANGED);
	Process(src);
	return 0;
}

int Process(IplImage *src)
{
	double *data_r, *u_r, *v_r, *c_r;
	double *data_g, *u_g, *v_g, *c_g;
	double *data_b, *u_b, *v_b, *c_b;
	int height, width, i, j, ka;
	double eps;
	int scale;
	IplImage *dst;
	
	dst = cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,3);

	eps=0.000001;
	height = src->height;
	width = src->width;
	
	//allocate memory for matrix
	data_r = (double *)malloc(sizeof(double)*height*width);
	u_r = (double *)malloc(sizeof(double)*height*height);
	v_r = (double *)malloc(sizeof(double)*width*width);

	data_g = (double *)malloc(sizeof(double)*height*width);
	u_g = (double *)malloc(sizeof(double)*height*height);
	v_g = (double *)malloc(sizeof(double)*width*width);

	data_b = (double *)malloc(sizeof(double)*height*width);
	u_b = (double *)malloc(sizeof(double)*height*height);
	v_b = (double *)malloc(sizeof(double)*width*width);

	memset(u_r,0,sizeof(double)*height*height);
	memset(v_r,0,sizeof(double)*width*width);

	memset(u_g,0,sizeof(double)*height*height);
	memset(v_g,0,sizeof(double)*width*width);

	memset(u_b,0,sizeof(double)*height*height);
	memset(v_b,0,sizeof(double)*width*width);

	//assign value
	for(i = 0;i < height;i++)
	{
		for(j = 0;j < width;j++) { data_r[i*width+j] = (double)(unsigned char)src->imageData[i*src->widthStep+j*src->nChannels+2];
			data_g[i*width+j] = (double)(unsigned char)src->imageData[i*src->widthStep+j*src->nChannels+1];
			data_b[i*width+j] = (double)(unsigned char)src->imageData[i*src->widthStep+j*src->nChannels+0];
		}
	}
	
	ka = max(height,width) + 1;
	bmuav(data_r, height, width, u_r, v_r, eps, ka);
	bmuav(data_g, height, width, u_g, v_g, eps, ka);
	bmuav(data_b, height, width, u_b, v_b, eps, ka);
	
	//dump svd, scale is selected by watching top xxx large data
	/*for (i=0; i<=100; i++)
	{ 
		for (j=0; j<=100; j++)
			printf("%f", data[i*width+j]);
		printf("\n");
	}*/
	
	//reconstruction
	scale = 50;
	for(i = scale;i<height;i++)
	{
		data_r[i*width+i] = 0;
		data_g[i*width+i] = 0;
		data_b[i*width+i] = 0;
	}

	/*c needs to be initilized here ,but in matrix mutiply funciton*/
	c_r = (double *)malloc(sizeof(double)*height*width);
	c_g = (double *)malloc(sizeof(double)*height*width);
	c_b = (double *)malloc(sizeof(double)*height*width);
	
	brmul(u_r, data_r ,height, height, width, c_r);
	brmul(c_r, v_r, height, width, width, data_r);

	brmul(u_g, data_g ,height, height, width, c_g);
	brmul(c_g, v_g, height, width, width, data_g);

	brmul(u_b, data_b ,height, height, width, c_b);
	brmul(c_b, v_b, height, width, width, data_b);

	//assign value
	for(i = 0;i < height;i++)
	{
		for(j = 0;j < width;j++) { dst->imageData[i*src->widthStep+j*src->nChannels + 0] = (unsigned char)data_b[i*width+j];
			dst->imageData[i*src->widthStep+j*src->nChannels + 1] = (unsigned char)data_g[i*width+j];
			dst->imageData[i*src->widthStep+j*src->nChannels + 2] = (unsigned char)data_r[i*width+j];
		}
	}
	cvSaveImage("result.jpg",dst);
	cvReleaseImage(&dst);
	free(u_r);free(v_r);free(c_r);free(data_r);
	free(u_g);free(v_g);free(c_g);free(data_g);
	free(u_b);free(v_b);free(c_b);free(data_b);
	return 0;
}

效果图:

  

注意:SVD压缩只是为了存储更少的数据来表达原始图像,在重构图像时,奇异值矩阵仍旧是要和原始图像大小一样的,只不过大部分地方用0填充罢了。

2 .     基于SVD的协同过滤推荐系统

这个例子比较有趣,推荐系统已存在很多年了。大家在网购时,商家总会根据大家的购买的历史记录给大家推荐新的商品及服务。实现方法有多种,这次只讲 基于协同过滤的方法。所谓协同过滤,就是将用户和其他用户的数据进行对比来实现推荐的。前端输入的原始数据如下:【评分:0-5】

【后话:很严重的问题,就是当用户间没有交集,推荐系统就无法工作咯】

    下面就构建一个餐饮网站的推荐系统。假设一个人在家决定外出吃饭,但是他并不知道该到哪去吃饭,该点什么菜。我们可以构建一个基本的推荐引擎,它能够帮助用户寻找没有尝过的菜肴,然后通过SVD来减少特征空间并提高反馈的速度。

a.      基本的推荐引擎

推荐系统的工作过程:给定一个用户,系统会为此用户返回N个最好的推荐,为了实现这一点我们需要做:

  •  寻找用户没有评级的菜肴
  •  在用户没有评级的物品中,对每个物品预计一个可能的评级分数。也就是说,我们的系统认为用户可能会对物品的打分
  •  对这些物品的评分从高到低排序,返回前N个

之前说将用户和其他用户的数据进行对比来实现推荐,那么我们借助什么手段预测评分呢?有两种方案:

基于用户的相似度:行与行之间的比较

基于菜肴的相似度:列与列之间的比较

【距离公式可采用欧氏距离、余弦距离、皮尔逊距离等等】

     那么到底使用哪一种相似度呢?这取决于菜肴和用户的数目。无论基于哪种相似度,它的计算时间都会随着用户/物品的数量的增加而增加。对于大部分产品导向的推荐引擎而言,用户的数量往往大于商品的数量。因此这里采用基于菜肴的相似度计算。

各种相似度计算代码:

void matrix_reverse(double *src,double *dest,int row,int col)
{
	int i,j;

	for(i = 0;i < col;i++)
	{
		for(j = 0;j < row;j++)
		{
			dest[i * row + j] = src[j * col + i];
		}
	}
}

/*欧氏距离*/
double ecludSim(double *dataMat,int *overLap, int n, int count, int item, int j)
{
	double total = 0.0;
	int i=0, index=0;

	for(i = 0;i < count;i++)
	{
		index = overLap[i];
		total = total + (dataMat[index*n+item] - dataMat[index*n+j]) * (dataMat[index*n+item] - dataMat[index*n+j]);
	}
	total = sqrt(total);

	return 1.0/(1.0 + total);
}
double ecludSim2(double *dataMat,int *overLap, int n, int scale, int item, int j)
{
	double total = 0.0;
	int i = 0;

	for(i = 0;i < scale;i++)
	{
		total = total + (dataMat[item*n+i] - dataMat[j*n+i]) * (dataMat[item*n+i] - dataMat[j*n+i]);
	}
	total = sqrt(total);

	return 1.0/(1.0 + total);
}

/*余弦距离*/
double cosSim(double *dataMat,int *overLap, int n, int count, int item, int j)
{
	double totalA=0.0, totalB=0.0, totalM=0.0;
	double result=0.0;
	int i, index;

	for(i = 0;i < count;i++)
	{
		index = overLap[i];
		totalA = totalA + dataMat[index*n+item] * dataMat[index*n+item];
		totalB = totalB + dataMat[index*n+j] * dataMat[index*n+j];
		totalM = totalM + dataMat[index*n+item] * dataMat[index*n+j];
	}
	result = totalM / (sqrt(totalA) * sqrt(totalB));

	return 0.5 + 0.5 * result;
}

double cosSim2(double *dataMat,int *overLap, int n, int scale, int item, int j)
{
	double totalA=0.0, totalB=0.0, totalM=0.0;
	double result=0.0;
	int i;

	for(i=0;i<scale;i++)
	{
		totalA = totalA + dataMat[item*n+i] * dataMat[item*n+i];
		totalB = totalB + dataMat[j*n+i] * dataMat[j*n+i];
		totalM = totalM + dataMat[item*n+i] * dataMat[j*n+i];
	}
	result = totalM / (sqrt(totalA) * sqrt(totalB));
	return 0.5 + 0.5 * result;
}

机器学习是用Python写的,操作矩阵很方便,这边C语言很笨重,所以没有做到代码重用,XXX2表示SVD推荐系统时用的距离函数,XXX是标准推荐系统】

预测评价的原理:假设要预测用户A的第k个未评价菜肴的评分,我们可以遍历整个数据矩阵寻找用户A和其他用户都评价过的其他菜肴j,计算用户A和其他用户针对菜肴j的评价的相似距离,若有多个这种共同的评价,则累加相似度。最后的计算公式:

第k个菜肴的评分 =  sum(累加的相似度 * A对j的评分 ) / 累加相似度

下面先给出主函数吧,否则太乱了:\

// recommand.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include 
#include 
#include 
#include "bmuav.c"
#include "brmul.c"
#define max(a,b)            (((a) > (b)) ? (a) : (b))

typedef double (*func)(double *dataMat,int *overLap, int n, int count, int item, int j);
typedef double (*Est)(double *dataMat, int user, func simMeas, int item, int m, int n);

extern void brmul(double *a,double *b,int m,int n,int k,double *c);
extern int bmuav(double *a,int m,int n,double *u,double *v,double eps,int ka);
extern void matrix_reverse(double *src,double *dest,int row,int col);

double standEst(double *dataMat, int user, func simMeas, int item, int m, int n);
double svdEst(double *dataMat, int user, func simMeas, int item, int m, int n);
double cosSim(double *dataMat,int *overLap, int n, int count, int item, int j);
double cosSim2(int *dataMat,int *overLap, int n, int scale, int item, int j);
double ecludSim(double *dataMat,int *overLap, int n, int count, int item, int j);
double ecludSim2(double *dataMat,int *overLap, int n, int count, int item, int j);
double recommend(double *dataMat,int n, int m, int user, func simMeas=cosSim, Est estMethod=standEst);

int _tmain(int argc, _TCHAR* argv[])
{
	int user;
	double data[7][5] = {{4, 4, 0, 2, 2},
					{4, 0, 0, 3, 3},
					{4, 0, 0, 1, 1},
					{1, 1, 1, 2, 0},
					{2, 2, 2, 0, 0},
					{5, 5, 5, 0, 0},
					{1, 1, 1, 0, 0}};
	double data_2[11][11] = {{2, 0, 0, 4, 4, 0, 0, 0, 0, 0, 0},
						{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5},
						{0, 0, 0, 0, 4, 0, 0, 1, 0, 4, 0},
						{3, 3, 4, 0, 3, 0, 0, 2, 2, 0, 0},
						{5, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0},
						{0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 0},
						{4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 5},
						{0, 0, 0, 0, 0, 4, 0, 0, 0, 0, 4},
						{0, 0, 0, 0, 0, 0, 5, 0, 0, 5, 0},
						{0, 0, 0, 3, 0, 0, 0, 0, 4, 5, 0},
						{1, 1, 2, 1, 1, 2, 1, 0, 4, 5, 0}};
	
	
	user = 2;

    //recommend(&data[0][0], 5, 7, user, cosSim, standEst);
	//recommend(&data[0][0], 5, 7, user, ecludSim, standEst);
	//recommend(&data_2[0][0], 11, 11, user, ecludSim2, svdEst);
	recommend(&data_2[0][0], 11, 11, user, cosSim2, svdEst);
	return 0;
}
double recommend(double *dataMat,int n, int m, int user, func simMeas, Est estMethod)
{
	int i,count, item,j;
	int *record=0;
	double *temp_vote=0;
	double temp=0, temp2=0;

	record = (int *)malloc(sizeof(int)*n);
	memset(record, 0, sizeof(int)*n);
	count = 0;

	//寻找user用户未评价物品
	for(i = 0;i < n;i++)
	{
		if(dataMat[user*n+i] == 0)
		{
			record[count++] = i;
		}
	}
	if (count == 0)
	{
		printf("该用户评价了所有的物品\n");
		return -1; //用户评价了所有的物品
	}

	temp_vote = (double *)malloc(sizeof(double)*count);
	memset(temp_vote,0,sizeof(double)*count);

	for(i=0;i<count;i++)
	{
		item = record[i];
		temp_vote[i] = estMethod(dataMat, user, simMeas, item, m, n);
	}

	//排序
	for(i=0;i<count;i++)
	{
		for(j=0;j

标准推荐系统:

double standEst(double *dataMat, int user, func simMeas, int item, int m, int n)
{
	double simTotal=0, ratSimTotal=0, userRating=0, similarity=0;
	int j=0, k=0, count=0;
	int *overLap=0;	//记录交集
	
	overLap = (int *)malloc(sizeof(int)*m);
	
	for (j = 0;j < n;j++)
	{
		userRating = dataMat[user*n+j];		//user用户评价过的
		if(userRating == 0)
			continue;
		
		count = 0;
		memset(overLap, 0, sizeof(int)*m);

		for(k = 0;k < m;k++) { if(dataMat[k*n+item] > 0 && dataMat[k*n+j] > 0)		//寻找用户都评级的两个物品
			{	
				overLap[count++] = k;
			}
		}

		if (count == 0)
		{
			similarity = 0;
		}else
		{
			similarity = simMeas(dataMat, overLap, n, count, item, j);
		}
		
		simTotal += similarity;
		ratSimTotal += similarity * userRating;
	}

	free(overLap);
	if(0 == simTotal)
	{
		return 0.0;
	}else
	{
		return ratSimTotal/simTotal;
	}
}

对于data矩阵,用户2,对应第三行,预测结果:

b.      基于SVD的推荐引擎

在数据矩阵非常稀疏时,基于SVD的推荐引擎性能就会比标准的好很多。我们利用SVD将所有菜肴隐射到一个低维空间中,再利用和前面一样的相似度计算方法来进行推荐。

double svdEst(double *dataMat, int user, func simMeas, int item, int m, int n)
{
	double *u, *v, *data_new, *I, *dataMat2, *dataMatCopy, *svdMat;
	double simTotal=0, ratSimTotal=0, userRating=0, similarity=0, eps=0;
	int i=0, j=0, k=0, count=0, ka=0,scale=0;

	u = (double *)malloc(sizeof(double)*m*m);
	v = (double *)malloc(sizeof(double)*n*n);
	dataMat2 = (double *)malloc(sizeof(double)*m*n);
	dataMatCopy = (double *)malloc(sizeof(double)*m*n);
	for(i=0;i<m*n;i++) 
		dataMatCopy[i] = dataMat[i];
	for(i=0;i<m*m;i++) 
		u[i] = 0;
	for(i=0;i<n*n;i++) 
		v[i] = 0;

	eps = 0.000001;
	ka = max(m,n) + 1;

	//奇异值分解
	bmuav(&dataMatCopy[0], m, n, u, v, eps, ka);

	//挑选合适的奇异值:打印出所有的,再挑选,此处打印略
	scale = 4;
	I = (double *)malloc(sizeof(double)*scale*scale);
	for(i=0;i<scale*scale;i++) 
		I[i] = 0.0;

	for(i = 0;i < scale;i++)
	{
		I[i*scale+i] = dataMatCopy[i*n+i];
		//printf("%f ", I[i*scale+i]);
	}
	//printf("\n");

	//将物品转换到低维空间,data_new = dataMat' * U[:,0:scale] * I
	data_new = (double *)malloc(sizeof(double)*n*scale);
	svdMat = (double *)malloc(sizeof(double)*n*scale);
	
	matrix_reverse(dataMat, dataMat2, m, n);
    brmul(dataMat2, u, n, m, scale, data_new);
	brmul(data_new, I, n, scale, scale, svdMat);

	for (j = 0;j < n;j++)
	{
		userRating = dataMat[user*n+j];		//user用户评价过的
		if((userRating == 0) || (j == item))
			continue;

		similarity = simMeas(svdMat, NULL, scale, scale, item, j);  //由于是SVD压缩的,不用考虑交集了		
		simTotal += similarity;
		ratSimTotal += similarity * userRating;
	}

	free(u);free(v);
	free(data_new);free(dataMat2);
	free(dataMatCopy);free(I);free(svdMat);

	if(0 == simTotal)
	{
		return 0.0;
	}else
	{
		return ratSimTotal/simTotal;
	}
	
}

对于矩阵 data_2,用户2,对应第三行,预测结果:

c.      现实中的挑战

(1)如何对推荐引擎进行评价

我们既没有预测的目标值,也没有用户来调查对他们预测结果的满意度。我们这时可以将已知的评价结果去掉,然后对他们进行预测,最后计算预测值和真实值之间 的差异。通常用于推荐引擎评价的指标是最小均方根误差,它 首先计算均方误差的平均值,然后取其平方根 【sqrt((sum((a-b)*(a-b)))/num)】

(2)真实的系统

这边代码为了保证可读性,没有效率:

a.我们不必在每次预测评分时都对数据矩阵进行SVD分解。在大规模数据集上,频繁进行SVD分解,只会拖慢效率。在真实系统中,SVD每天运行一次,并且还要离线运行。

b.规模扩展性的挑战:

数据矩阵很稀疏,我们有很多0,我们能否只存储非零元素来节省内存计算开销。

相似度计算也导致了计算资源浪费,每次需要一个推荐得分时,都要计算很多物品的相似度得分,这些相似度得分能否被其他用户重复使用?实际系统中,会进行离线计算,并保存相似度得分。

c.如何在缺乏数据时,给出好的推荐?

实际系统将推荐系统看成是搜索问题,我们可能要使用需要推荐物品的属性。在上述餐馆例子里,我们可以通过各种标签来标记菜肴,比如素食、美式烤肉、价格贵等。同时我们也可以将这些属性作为相似度计算所需要的数据。这个被称为基于内容的推荐。

所有代码下载地址:

温馨提示:此处内容需要评论本文并通过审核后,刷新本页才能查看,请勿使用垃圾评论!
如果查看不了,删除电脑缓存,重新刷新本页
游客回复后依然不可见!

CV视觉网,转载请注明来源:CV视觉网 » 机器学习—奇异值分解

分享到:更多 ()
微信 OR 支付宝 扫描二维码
为本文作者 打个赏
pay_weixinpay_weixin
金额随意 快来“打”我呀~

请选择你看完该文章的感受:

0不错 1超赞 0无聊 0扯淡 0不解 0路过