Singular Value Decomposition

詹令
lealzhan@126.com
2016.1.17

SVD是什么

TODO
x¹²

SVD的数值求解

TODO

numerical receipe³:

#define	MIN(a,b)		((a)<(b)?(a):(b))
#define	MAX(a,b)		((a)>(b)?(a):(b))
#define SIGN(a)			((a)<0?-1:1)

///
//  SVD function <from numerical recipes in C++>
//		Given a matrix a[1..m][1..n], this routine computes its singular value
//		decomposition, A = U.W.VT.  The matrix U replaces a on output.  The diagonal
//		matrix of singular values W is output as a vector w[1..n].  The matrix V (not
//		the transpose VT) is output as v[1..n][1..n].
///
 double pythag(double a, double b)
{
	double at = fabs(a), bt = fabs(b), ct, result;
	if (at > bt)       { ct = bt / at; result = at * sqrt(1.0 + ct * ct); }
	else if (bt > 0.0) { ct = at / bt; result = bt * sqrt(1.0 + ct * ct); }
	else result = 0.0;
	return(result);
}	

 void SVD(double* u, int m, int n, double* w, double* v)
{
	bool flag;
	int i,its,j,jj,k,l,nm;
	double anorm,c,f,g,h,s,scale,x,y,z;
	double *rv1=new double[n];
	g = scale = anorm = 0.0; //Householder reduction to bidiagonal form.
	for (i=0;i<n;i++) 
	{
		l=i+1;
		rv1[i]=scale*g;
		g=s=scale=0.0;
		if (i < m) 
		{
			for (k=i;k<m;k++) scale += fabs(u[k*n+i]);
			if (scale != 0.0) 
			{
				for (k=i;k<m;k++) 
				{
					u[k*n+i] /= scale;
					s += u[k*n+i]*u[k*n+i];
				}
				f=u[i*n+i];
				g = -sqrt(s)*SIGN(f);
				h=f*g-s;
				u[i*n+i]=f-g;
				for (j=l;j<n;j++) 
				{
					for (s=0.0,k=i;k<m;k++) s += u[k*n+i]*u[k*n+j];
					f=s/h;
					for (k=i;k<m;k++) u[k*n+j] += f*u[k*n+i];
				}
				for (k=i;k<m;k++) u[k*n+i] *= scale;
			}
		}
		w[i]=scale *g;
		g=s=scale=0.0;
		if (i+1 <= m && i+1 != n) 
		{
			for (k=l;k<n;k++) scale += fabs(u[i*n+k]);
			if (scale != 0.0) 
			{
				for (k=l;k<n;k++) 
				{
					u[i*n+k] /= scale;
					s += u[i*n+k]*u[i*n+k];
				}
				f=u[i*n+l];
				g = -sqrt(s)*SIGN(f);
				h=f*g-s;
				u[i*n+l]=f-g;
				for (k=l;k<n;k++) rv1[k]=u[i*n+k]/h;
				for (j=l;j<m;j++) 
				{
					for (s=0.0,k=l;k<n;k++) s += u[j*n+k]*u[i*n+k];
					for (k=l;k<n;k++) u[j*n+k] += s*rv1[k];
				}
				for (k=l;k<n;k++) u[i*n+k] *= scale;
			}
		}
		anorm=MAX(anorm,(fabs(w[i])+fabs(rv1[i])));
	}
	for (i=n-1;i>=0;i--) 
	{ //Accumulation of right-hand transformations.
		if (i < n-1) 
		{
			if (g != 0.0) 
			{
				for (j=l;j<n;j++) //Double division to avoid possible under
					v[j*n+i]=(u[i*n+j]/u[i*n+l])/g;
				for (j=l;j<n;j++) 
				{
					for (s=0.0,k=l;k<n;k++) s += u[i*n+k]*v[k*n+j];
					for (k=l;k<n;k++) v[k*n+j] += s*v[k*n+i];
				}
			}
			for (j=l;j<n;j++) v[i*n+j]=v[j*n+i]=0.0;
		}
		v[i*n+i]=1.0;
		g=rv1[i];
		l=i;
	}
	for (i=MIN(m,n)-1;i>=0;i--) 
	{	//Accumulation of left-hand transformations.
		l=i+1;
		g=w[i];
		for (j=l;j<n;j++) u[i*n+j]=0.0;
		if (g != 0.0) 
		{
			g=1.0/g;
			for (j=l;j<n;j++) 
			{
				for (s=0.0,k=l;k<m;k++) s += u[k*n+i]*u[k*n+j];
				f=(s/u[i*n+i])*g;
				for (k=i;k<m;k++) u[k*n+j] += f*u[k*n+i];
			}
			for (j=i;j<m;j++) u[j*n+i] *= g;
		} 
		else for (j=i;j<m;j++) u[j*n+i]=0.0;
		++u[i*n+i];
	}
	for (k=n-1;k>=0;k--) 
	{ //Diagonalization of the bidiagonal form: Loop over
		for (its=0;its<30;its++) 
		{ //singular values, and over allowed iterations.
			flag=true;
			for (l=k;l>=0;l--) 
			{ //Test for splitting.
				nm=l-1;
				if ((double)(fabs(rv1[l])+anorm) == anorm) 
				{
					flag=false;
					break;
				}
				if ((double)(fabs(w[nm])+anorm) == anorm) break;
			}
			if (flag) 
			{
				c=0.0; //Cancellation of rv1[l], if l > 0.
				s=1.0;
				for (i=l;i<k+1;i++) 
				{
					f=s*rv1[i];
					rv1[i]=c*rv1[i];
					if ((double)(fabs(f)+anorm) == anorm) break;
					g=w[i];
					h=pythag(f,g);
					w[i]=h;
					h=1.0/h;
					c=g*h;
					s = -f*h;
					for (j=0;j<m;j++) 
					{
						y=u[j*n+nm];
						z=u[j*n+i];
						u[j*n+nm]=y*c+z*s;
						u[j*n+i]=z*c-y*s;
					}
				}
			}
			z=w[k];
			if (l == k) 
			{	// Convergence.
				if (z < 0.0) 
				{ //Singular value is made nonnegative.
					w[k] = -z;
					for (j=0;j<n;j++) v[j*n+k] = -v[j*n+k];
				}
				break;
			}
			if (its == 29) printf("Error: no convergence in 30 svdcmp iterations");
			x=w[l]; //Shift from bottom 2-by-2 minor.
			nm=k-1;
			y=w[nm];
			g=rv1[nm];
			h=rv1[k];
			f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
			g=pythag(f,1.0);
			f=((x-z)*(x+z)+h*((y/(f+fabs(g)*SIGN(f)))-h))/x;
			c=s=1.0; //Next QR transformation:
			for (j=l;j<=nm;j++) 
			{
				i=j+1;
				g=rv1[i];
				y=w[i];
				h=s*g;
				g=c*g;
				z=pythag(f,h);
				rv1[j]=z;
				c=f/z;
				s=h/z;
				f=x*c+g*s;
				g=g*c-x*s;
				h=y*s;
				y *= c;
				for (jj=0;jj<n;jj++) 
				{
					x=v[jj*n+j];
					z=v[jj*n+i];
					v[jj*n+j]=x*c+z*s;
					v[jj*n+i]=z*c-x*s;
				}
				z=pythag(f,h);
				w[j]=z; //Rotation can be arbitrary if z D 0.
				if (z) 
				{
					z=1.0/z;
					c=f*z;
					s=h*z;
				}
				f=c*g+s*y;
				x=c*y-s*g;
				for (jj=0;jj<m;jj++) 
				{
					y=u[jj*n+j];
					z=u[jj*n+i];
					u[jj*n+j]=y*c+z*s;
					u[jj*n+i]=z*c-y*s;
				}
			}
			rv1[l]=0.0;
			rv1[k]=f;
			w[k]=x;
		}
	}
	delete []rv1;
}

python⁴ :
U, Sigma, VT = linalg.svd(data);

应用

求解线性方程组

比如在numerical receipe³中，用SVD求解线性方程组的代码如下：

///
//  SVD function <from numerical recipes in C++>
//		Given a matrix a[1..m][1..n], this routine computes its singular value
//		decomposition, A = U.W.VT.  The matrix U replaces a on output.  The diagonal
//		matrix of singular values W is output as a vector w[1..n].  The matrix V (not
//		the transpose VT) is output as v[1..n][1..n].
///
double pythag(double a, double b)

void SVD_decompose(double* u, int m, int n, double* w, double* v)

void SVD_reorder(double* u, int m, int n, double* w, double* v)

//Solve A  x D b for a vector x using the pseudoinverse of A as obtained by SVD.If positive,
//thresh is the threshold value below which singular values are considered as zero.If thresh is
//negative, a default based on expected roundoff error is used.
void SVD_solve(double* u, int m, int n, double* w, double* v, double* x, double* b, double thresh/* = -1.*/)

void SVD_linear_equation_main(double* A, int m, int n, double* X, double* B)
{
	double* u = A;
	double* w = (double*)malloc( sizeof(double) * n );
	double* v = (double*)malloc( sizeof(double) * n * n);

	SVD_decompose(u, m, n, w, v);
	SVD_reorder(u, m, n, w, v);
	SVD_solve(u, m, n, w, v, X, B, -1);
}

在实际应用中，我曾经在图像的透视变换中使用了SVD来求解其涉及到的低维度的线性方程组。

Latent Semantic Analysis 模型, 隐性语义分析

动力学模拟

在Least-Squares Rigid Motion Using SVD⁵ 中，求解刚体运动的问题转换成了对协方差矩阵进行SVD在得到变换矩阵的问题。具体见⁵ 以及 博文

Deformation Gradient Matrix 经过SVD（polar decomposiiton），得到R*U，R是旋转矩阵，U代表Scale,strain, crack等。

数据压缩 / 降维 / 除噪

原始矩阵svd分解后，如上图。仅需存储黑色部分矩阵。将黑色部分矩阵相乘后即可近似还原原始矩阵。一般情况下，采用能量90%原则决定近似程度（黑色矩阵大小)⁴。

机器学习实战⁴ 14.5 图像压缩
数学之美⁶ 数据分析。降维

张正友相机校正

基于协同过滤的推荐引擎

[5] 14.4
协同过滤 collaborative filtering ：通过将_用户_和_其他用户_的数据进行对比来实现推荐。
–待补充

思考小问题

SVD得到奇异值a（中间那个矩阵的对角线上）和被分解矩阵的特征值b（PCA）之间的关系：
a=sqrt(b)?
两者概念上有何不同？

参考

---

1. 格罗布. 矩阵计算[M]. 科学出版社, 2001. ↩︎

2. 奇异值分解（SVD）详解, CSDN BLOG LINK ↩︎

3. Numerical receipes in C: The art of scientific computing, 2nd Ed ↩︎ ↩︎

4. 机器学习实战 ↩︎ ↩︎ ↩︎

5. Olga Sorkine-Hornung, Least-Squares Rigid Motion Using SVD ↩︎ ↩︎

6. 吴军，数学之美, 矩阵分类和文本处理中的两个基本问题 ↩︎