/** C++ function for SVD
函数原型:
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是输入矩阵，假设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，0<K<=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=0i<ni++)
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=0i<ni++)
x[i]/=r
return r
}
inline double product(double*a, double *b,int n){
double r=0
for(int i=0i<ni++)
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=0i<ni++)
b[i]-=r*a[i]
}
bool svd(vector<vector<double>>A, int K, std::vector<std::vector<double>>&U, std::vector<double>&S, std::vector<std::vector<double>>&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=0i<Ki++)
U[i].resize(M,0)
V.resize(K)
for(int i=0i<Ki++)
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=0i<Mi++)
橡灶 left_vector[i]= (float)rand() / RAND_MAX
if(normalize(left_vector, M)>eps)
break
}
int col=0
for(int col=0col<Kcol++){
double diff=1
double r=-1
for(int iter=0diff>=eps &&iter<MAX_ITERiter++){
兄纳 memset(next_left_vector,0,sizeof(double)*M)
memset(next_right_vector,0,sizeof(double)*N)
for(int i=0i<Mi++)
for(int j=0j<Nj++)
next_right_vector[j]+=left_vector[i]*A[i][j]
r=normalize(next_right_vector,N)
if(r<eps) break
for(int i=0i<coli++)
orth(&V[i][0],next_right_vector,N)
normalize(next_right_vector,N)
for(int i=0i<Mi++)
for(int j=0j<Nj++)
next_left_vector[i]+=next_right_vector[j]*A[i][j]
r=normalize(next_left_vector,M)
if(r<eps) break
for(int i=0i<coli++)
orth(&U[i][0],next_left_vector,M)
normalize(next_left_vector,M)
diff=0
for(int i=0i<Mi++){
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(vector<vector<double>>&A){
for(int i=0i<A.size()i++){
for(int j=0j<A[i].size()j++){
cout<<setprecision(3)<<A[i][j]<<' '
}
cout<<endl
}
}
int main(){
int m=10
int n=5
srand(time(0))
vector<vector<double>>A
A.resize(m)
for(int i=0i<mi++){
A[i].resize(n)
for(int j=0j<nj++)
A[i][j]=(float)rand()/RAND_MAX
}
print(A)
cout<<endl
vector<vector<double>>U
vector<double>S
vector<vector<double>>V
svd(A,2,U,S,V)
cout<<"U="<<endl
print(U)
cout<<endl
cout<<"S="<<endl
for(int i=0i<S.size()i++){
cout<<S[i]<<' '
}
cout<<endl
cout<<"V="<<endl
print(V)
return 0
}
/*
本程序在linux g++下编译通过
bool svd(vector<vector<double>>A, int K, vector<vector<double>>&U, vector<double>&S, vector<vector<double>>&V)
A: 输入待分解矩阵
K: 输入，取前K大奇异值及奇异向量
U[0],U[1],...,U[K-1]: 前K大奇异值对应的左奇异向量
S[0],S[1],...,S[K-1]: 前K大奇异值 S[0]>=S[1]>=...>=S[K-1]
V[0],V[1],...,V[K-1]: 前K大奇异值对应的右奇异向量
*/
#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=0i<ni++)
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=0i<ni++)
x[i]/=r
return r
}
inline double product(double*a, double *b,int n){
double r=0
for(int i=0i<ni++)
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=0i<ni++)
b[i]-=r*a[i]
}
bool svd(vector<vector<double>>A, int K, vector<vector<double>>&U, vector<double>&S, vector<vector<double>>&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=0i<Ki++)
U[i].resize(M,0)
V.resize(K)
for(int i=0i<Ki++)
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]
int col=0
for(int col=0col<Kcol++){
double diff=1
double r=-1
while(1){
for(int i=0i<Mi++)
left_vector[i]= (float)rand() / RAND_MAX
if(normalize(left_vector, M)>eps)
break
}
for(int iter=0diff>=eps &&iter<MAX_ITERiter++){
memset(next_left_vector,0,sizeof(double)*M)
memset(next_right_vector,0,sizeof(double)*N)
for(int i=0i<Mi++)
for(int j=0j<Nj++)
next_right_vector[j]+=left_vector[i]*A[i][j]
r=normalize(next_right_vector,N)
if(r<eps) break
for(int i=0i<coli++)
orth(&V[i][0],next_right_vector,N)
normalize(next_right_vector,N)
for(int i=0i<Mi++)
for(int j=0j<Nj++)
next_left_vector[i]+=next_right_vector[j]*A[i][j]
r=normalize(next_left_vector,M)
if(r<eps) break
for(int i=0i<coli++)
orth(&U[i][0],next_left_vector,M)
normalize(next_left_vector,M)
diff=0
for(int i=0i<Mi++){
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{
cout<<r<<endl
break
}
}
delete [] next_left_vector
delete [] next_right_vector
delete [] left_vector
delete [] right_vector
return true
}
void print(vector<vector<double>>&A){
}
int main(){
int m=10
int n=8
int k=5
//分解一个10*8的矩阵A，求其前5个奇异值和奇异向量
srand(time(0))
vector<vector<double>>A
A.resize(m)
for(int i=0i<mi++){
A[i].resize(n)
for(int j=0j<nj++)
A[i][j]=(float)rand()/RAND_MAX-0.5
}
cout<<"A="<<endl
for(int i=0i<A.size()i++){
for(int j=0j<A[i].size()j++){
cout<<setw(12)<<A[i][j]<<' '
}
cout<<endl
}
cout<<endl
vector<vector<double>>U
vector<double>S
vector<vector<double>>V
svd(A,k,U,S,V)
cout<<"U="<<endl
for(int i=0i<U[0].size()i++){
for(int j=0j<U.size()j++){
cout<<setw(12)<<U[j][i]<<' '
}
cout<<endl
}
cout<<endl
cout<<"S="<<endl
for(int i=0i<S.size()i++){
cout<<setw(7)<<S[i]<<' '
}
cout<<endl
cout<<"V="<<endl
for(int i=0i<V[0].size()i++){
for(int j=0j<V.size()j++){
cout<<setw(12)<<V[j][i]<<' '
}
cout<<endl
}
return 0
}
211 SVD算法SVD算法可用来求解大多数的线性最小二乘法问题. SVD 算法基于如下分解定理:对任
意的矩阵 Am ×n ,当其行数 m 大于等于列数 n 时,可以分解为正交矩阵 Um ×n , 非负对角矩阵卖枯
Wn×n以及正交矩阵Vn×n的转置的乘积,即
Am×n = Um×n ·[diag( wj
) ] n×n ·V T
n×n , (2)
其中 wj ≥ 0 (1 ≤j ≤n) U , V 为正交矩阵,即满足
6
m
i =1
uijuik = δ jk ,
6
n
i =1
vijvik = δ jk
(1 ≤j , k ≤n) , (3)当 m <n 时,SVD也可以执行,在这种情况下,奇异值 wj = 0 ( m + 1 ≤j ≤n) ,并且 U 中相应的
列都是零,这时(3)式仅对 j , k ≤m 时成立. 故不管矩阵 A 是否是奇异, (2)式的分解总可以进
行,而且这个分解几乎是惟一的. 也就是说,其分解形式惟一到:对矩阵 U 的列、 W 的元素和
V 的列能做相同的置换,或者矩阵 U 和V 的任意列的线性组合,在 W 中对应的元素仍恰好完
全相同.
SVD分解明中渣洞确地构造了矩阵零空间和值域的正交标准化基. 特别地,对 U 的列,若与其
标号相同的元素 wj 为非零元,则其列为值域的一个正交标准化梁迹的基础矢量对 V 的列,若与
其标号相同的 wj 为零,则其列为零空间的一个正交标准化基. 对于如下的多指数衰减 T2 模
型, 有
y = M ·f , (4)
其中 y = ( y1 , y2 , …, yn )
T
为测量的自旋回波衰减信号, M = [ mij ] n ×m = [ exp ( - ti/
T2 j
) ] n ×m f = ( f 1 , f 2 , …, f m)
T
为弛豫时间 T2 j对应的各点的幅度值, T2 j
( j = 1 ,2 , …, m)为预先
指定的 T2 时间分布系列,典型的取法为在( T2min , T2max)区间内对数均匀地选取 m 个点,我们
称为弛豫时间布点,也可采用2的幂指数布点、 线性均匀布点等方式. 矩阵条件数的定义为矩
阵的最大特征值与最小特征值的比值. 若矩阵的条件数为无穷大,则该矩阵奇异若矩阵的条
件数太大,即其倒数超出了机器的浮点精度,则称该矩阵为病态的矩阵. 若直接采用 G auss分
解求上式,几乎是不可能的,原因是矩阵 M 的条件数相当大,例如:若回波间隔τ= 1. 2 ×10 - 3
s , T2 在0. 1 × 10 - 3
～10000 ×10 - 3
s内对数均匀地取50 个点,则矩阵 M 的条件数可达1016
数
量级,很明显,矩阵 M 是高度病态的. 采用 SVD 分解法来求解上式,系数矩阵 Mn×m =
Un×m ·[diag( wj
) ] m×m · VT
m×m , 这里 U , V 为正交矩阵,diag( wj
)为对角矩阵,其对角元递减排
列,则我们就可以很容易地求得最小二乘意义下的解为
^ f = V · diag
1
w1
,
1
w2
, …,
SNR
w1
, 0 , …,0 ·( UT
·y) . (5)
这里给出了矩阵条件数小于等于SNR的限制,避免了解的不稳定性. 其中 SNR为从测量数据
中估计出的信噪比. SNR定义为第1个回波的幅度值除以误差矢量 r ( r = y - M· ^ f )的标准差
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）