基于单边jacobi的奇异值分解（SVD）

基于单边jacobi的奇异值分解（SVD）

对于奇异值分解（SVD），相信很多学过线性代数/高等代数的同学都不会很陌生，但是怎么实现呢？接下来就来详细说说。

　　为了方便讨论，本文所有的讨论仅限定于实数空间。奇异值分解的含义就是将矩阵A分解成一个酉矩阵U，一个准对角矩阵S还有一个酉矩阵V。首先要说明一点，SVD分解是存在的但不唯一，这个有兴趣的读者可以思考一下。利用数学公式可以写成：

A=USV′ A∈M∗N U∈M∗M S∈M∗N V∈N∗N\ A= USV'\quad\ A\in\mathbb M*N\quad\ U\in\mathbb M*M\quad\ S\in\mathbb M*N\quad\ V\in\mathbb N*N\quad

　　其实大家关心的是，如何将一个实矩阵分解成这样的形式，以及为什么要这样做。了解PCA（主成分分析）的同学肯定知道，PCA就是为了降维的，那么SVD也是做这样的事情的，只要求出准对角阵S，就可以把无关紧要的几行或者几列去掉了。这在实际应用中是有很大的价值的。

　　

　　了解完为什么要这样做以后，就来看看实现了。为了后面并行的文章，这里不讨论QR分解或者随机投影等分解方法。基于jacobi进行SVD主要是有两种方法，如下表所示，一种是双边jacobi，另一种是单边jacobi。

方法	实现平台
双边jacobi	MATLAB
单边jacobi	MATLAB

1、双边jacobi

　　其实这个思想很巧妙也很简单（这就是数学的魅力），由于一个矩阵左乘或者右乘一个单位阵仅仅是对行或者列做了变换，恰巧这些单位阵又都是正交的，正交阵与正交阵相乘还是正交的。于是Ｕ和Ｖ就可以被求出来了。多说无益，贴几张图片就可以了。

　 　　　　　


　 　　　　　

　 　　为了把Ｂ的非对角元素化为０

　 　　　　　


　 　　于是就可以求出相应的ｃ和ｓ了。

　 　　


　 　　

2、单边jacobi

　　单边jacobi的思想更加巧妙，事实上由于 A=USV′\ A= USV'，则有 AV=US\ AV= US，我们观察到等式右边是一个正交阵，因此一个基本的思想就是将 A\ A化成正交阵。

　　


　　

3、实现

　　为了方便快捷，实现的平台选择MATLAB，当然读者也可以选择C++或者其他语言去做。由于每次清扫的过程中都涉及到旋转变换，因此要保证每次清扫所有列都能被清扫并且不重复。于是，我们可以采用round_robin或者前向序列（后向序列）去做这个事情。

　　round_robin序列

　　


　　ring序列

　　

　　　


　　

　　根据以上描述，我们可以得到单边jacobi的伪代码：

　　


　　需要进一步解释的是伪代码中的P指的就是把 A\ A的列范数从大到小排列后所得的index。为什么每次都要将 A\ A的列范数从大到小排列呢？这是因为，旋转以后列范数大的那一列列范数还会增大，而列范数小的那一列就会减小，所以为了保证列范数不会震荡，这里面就有一个小技巧了。我们就保证每次进行变换的两列总是序号在前（因为序号在前的列范数大）的那一列增大。但是这里又有一个问题： Aro∗Aco\ Aro*Aco这个数不总是大于0的，关于这一点，希望读者自己思考了。不过其实最后的实验证明这一步做或者不做都不太重要。

　　

function [ A,V ] = test_hestenes_jacobi( A )
%主程序部分
[n,m]=size(A);
B=zeros(n,m);
n_element=n*(n-1)/2;
n_sweep=0;
e=1e-8;
convergence=0;
d=zeros(n,1);
V=eye(n);
for i=1:n
d(i)=A(1:n,i)'*A(1:n,i);
end
% [d_sort,P]=sort(d,'descend');
% for i=1:m
%     B(:,i)=A(:,P(i));
% end
% A=B;
% for i=1:n
%     d(i)=A(1:n,i)'*A(1:n,i);
% end
for i=1:30
n_sweep=n_sweep+1;
n_clear=0;
[d_sort,P]=sort(d,'descend');
P_a=P(n-1:-2:1);
P_b=P(n:-2:2);
count_i=1;
for k=1:n-1
for j=1:n/2
ro=P_a(j);
co=P_b(j);
if(ro>co)  %为了保证序号小列的范数可能增大（可以不做）
temp=ro;
ro=co;
co=temp;
end
Aro=A(:,ro);
Aco=A(:,co);
d_norm=Aro'*Aco;
error=n*e*sqrt(d(ro)*d(co));
if(abs(d_norm)>error)
if(d(ro)<d(co)) %如果序号小的列范数小，那么就交换两列（可以不做）。
exchange=1;
%                   temp=Aro;
%                   Aro=Aco;
%                   Aco=temp;

temp=d(ro);
d(ro)=d(co);
d(co)=temp;

%                   temp=V(:,ro);
%                   V(:,ro)=V(:,co);
%                   V(:,co)=temp;
else
exchange=0;
end
if(exchange==0)
[Aro,Aco,dro,dco,Vro,Vco]=one_side_rotate(Aro,Aco,d(ro),d(co),d_norm,V(:,ro),V(:,co));%进行旋转
else
[Aro,Aco,dro,dco,Vro,Vco]=one_side_rotate(Aco,Aro,d(ro),d(co),d_norm,V(:,co),V(:,ro)); %进行旋转
end
A(:,ro)=Aro;
A(:,co)=Aco;
d(ro)=dro;
d(co)=dco;
V(:,ro)=Vro;
V(:,co)=Vco;
else
n_clear=n_clear+1;
end
if(n_clear==n_element)
convergence=1;
break;
end
end
if(convergence==1)
break;
end
[P_a,P_b,count_i] = initilize_index_of_one_side_jacobi(P_a,P_b,count_i,k);%利用ring序列前向序列调整序列号，保证每列都能旋转。
end
if(convergence==1)
break;
end
end
n_sweep
end

function [ A1,A2,d1,d2,V1,V2 ] = one_side_rotate(Aro,Aco,dro,dco,d,Vro,Vco)
%旋转变换
ct=(dro-dco)/(2*d);
t=sign(ct)/(abs(ct)+sqrt(1+ct^2));
c=1/sqrt(1+t^2);
s=t*c;
A1=c*Aro+s*Aco;
A2=-s*Aro+c*Aco;
d1=A1'*A1;
d2=A2'*A2;%在C++里面这一步非常耗时间，可以利用d1=c^2*d1+s^2*d2+2c*s*d;
V1=c*Vro+s*Vco;
V2=-s*Vro+c*Vco;
end

function [ P_a,P_b,count_i] = initilize_index_of_one_side_jacobi(P_a,P_b,i,k)
%这一步主要是做ring序列的变换
if(k==(2*i+1))
i=i+1;
end
[P_a,P_b]=push_forward(P_a,P_b,i);
count_i=i;
end

function [a,b] = push_forward( a,b,i )
%UNTITLED6 Summary of this function goes here
%   Detailed explanation goes here
temp=a(i);
a(i)=b(i);
b(i)=temp;
n=length(b);
temp=b(n);
for i=n:-1:2
b(i)=b(i-1);
end
b(1)=temp;

end

至此，整个单边JACOBI的程序就给出来了，有兴趣的读者可以把ring序列换成round_robin测试一下到底那个序列的效果比较好。下面给出round_robin序列的代码。

function [ a,b ] = round_robin(a,b)
%generate round robin sequence
% update P_a and P_b
n=length(a);
temp_bn=b(n);
temp_a_1=a(1);
for i=n:-1:2
b(i)=b(i-1);
end

for i=1:n-2
a(i)=a(i+1);
end
a(n-1)=temp_bn;
b(1)=temp_a_1;

end

[1]: 唐吉卓. 基于GPU平台的SVD并行计算研究与实现[D]. 电子科技大学, 2014.