/** C++ function for SVD
站在用户的角度思考问题,与客户深入沟通,找到如皋网站设计与如皋网站推广的解决方案,凭借多年的经验,让设计与互联网技术结合,创造个性化、用户体验好的作品,建站类型包括:成都网站建设、网站制作、企业官网、英文网站、手机端网站、网站推广、空间域名、虚拟主机、企业邮箱。业务覆盖如皋地区。
函数原型:
bool svd(vectorvectordouble A, int K, std::vectorstd::vectordouble U, std::vectordouble S, std::vectorstd::vectordouble V);
其中
A是输入矩阵,假设A的维数是m*n,那么本函数将A分解为U diag(S) V'
其中U是m*K的列正交的矩阵. V是n*K的列正交矩阵,S是K维向量。K由第二个参数指定。
U的第i列是A的第i大奇异值对应的左歧义向量,S[i]=A的第 i大奇异值,V的第i列是A的第i大奇异值对应的右歧义响亮.
K是需要分解的rank,0K=min(m,n)
本程序采用的是最基本幂迭代算法,在linux g++下编译通过
**/
#include cmath
#include iostream
#include iomanip
#include cstdlib
#include cstring
#include fstream
#include vector
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;in;i++)
r+=x[i]*x[i];
return sqrt(r);
}
double normalize(double *x, int n){
double r=get_norm(x,n);
if(reps)
return 0;
for(int i=0;in;i++)
x[i]/=r;
return r;
}
inline double product(double*a, double *b,int n){
double r=0;
for(int i=0;in;i++)
r+=a[i]*b[i];
return r;
}
void orth(double *a, double *b, int n){//|a|=1
double r=product(a,b,n);
for(int i=0;in;i++)
b[i]-=r*a[i];
}
bool svd(vectorvectordouble A, int K, std::vectorstd::vectordouble U, std::vectordouble S, std::vectorstd::vectordouble V){
int M=A.size();
int N=A[0].size();
U.clear();
V.clear();
S.clear();
S.resize(K,0);
U.resize(K);
for(int i=0;iK;i++)
U[i].resize(M,0);
V.resize(K);
for(int i=0;iK;i++)
V[i].resize(N,0);
srand(time(0));
double *left_vector=new double[M];
double *next_left_vector=new double[M];
double *right_vector=new double[N];
double *next_right_vector=new double[N];
while(1){
for(int i=0;iM;i++)
left_vector[i]= (float)rand() / RAND_MAX;
if(normalize(left_vector, M)eps)
break;
}
int col=0;
for(int col=0;colK;col++){
double diff=1;
double r=-1;
for(int iter=0;diff=eps iterMAX_ITER;iter++){
memset(next_left_vector,0,sizeof(double)*M);
memset(next_right_vector,0,sizeof(double)*N);
for(int i=0;iM;i++)
for(int j=0;jN;j++)
next_right_vector[j]+=left_vector[i]*A[i][j];
r=normalize(next_right_vector,N);
if(reps) break;
for(int i=0;icol;i++)
orth(V[i][0],next_right_vector,N);
normalize(next_right_vector,N);
for(int i=0;iM;i++)
for(int j=0;jN;j++)
next_left_vector[i]+=next_right_vector[j]*A[i][j];
r=normalize(next_left_vector,M);
if(reps) break;
for(int i=0;icol;i++)
orth(U[i][0],next_left_vector,M);
normalize(next_left_vector,M);
diff=0;
for(int i=0;iM;i++){
double d=next_left_vector[i]-left_vector[i];
diff+=d*d;
}
memcpy(left_vector,next_left_vector,sizeof(double)*M);
memcpy(right_vector,next_right_vector,sizeof(double)*N);
}
if(r=eps){
S[col]=r;
memcpy((char *)U[col][0],left_vector,sizeof(double)*M);
memcpy((char *)V[col][0],right_vector,sizeof(double)*N);
}else
break;
}
delete [] next_left_vector;
delete [] next_right_vector;
delete [] left_vector;
delete [] right_vector;
return true;
}
void print(vectorvectordouble A){
for(int i=0;iA.size();i++){
for(int j=0;jA[i].size();j++){
coutsetprecision(3)A[i][j]' ';
}
coutendl;
}
}
int main(){
int m=10;
int n=5;
srand(time(0));
vectorvectordouble A;
A.resize(m);
for(int i=0;im;i++){
A[i].resize(n);
for(int j=0;jn;j++)
A[i][j]=(float)rand()/RAND_MAX;
}
print(A);
coutendl;
vectorvectordouble U;
vectordouble S;
vectorvectordouble V;
svd(A,2,U,S,V);
cout"U="endl;
print(U);
coutendl;
cout"S="endl;
for(int i=0;iS.size();i++){
coutS[i]' ';
}
coutendl;
cout"V="endl;
print(V);
return 0;
}
数据挖掘比赛算法
examples/src/main/java/org/apache/mahout/cf/taste/example/kddcup/track1/svd
推荐系统中利用SVD实现降维
core/src/main/java/org/apache/mahout/cf/taste/impl/recommender/svd
String getStr(String str) {
if (str.length() = 2) {
return str;
}
return str.substring(0, 2) + "#" + getStr(str.substring(2));
}
不要太简单
SVD相当于:
C为mxn阶矩阵,U为mxm阶酉矩阵,E为mxn阶实数对角矩阵,V为nxn阶酉矩阵。
E矩阵对角线上的元素等于C的奇异值。
在OpenCV中可以用
solve(InputArray src1, InputArray src2, OutputArray dst, intflags=DECOMP_SVD)
你是不是要用SVD求解最小二乘问题?
线性方程组Cx=b,求其最小二乘解。
你可以先求出C的伪逆C+。
C+=V(E+)(UT)
(E+)是E的伪逆,将E主对角线上每个非零元素都求倒数之后再转置得到.
(UT)是U的转置。
最后计算出最小二乘解为:(C+)b