矩阵奇异值计算的一种新方法——基于R语言实现

Runnn1

已于 2022-10-10 17:07:49 修改

阅读量1k

点赞数 3

文章标签： r语言

于 2022-10-10 17:00:43 首次发布

本文链接：https://blog.csdn.net/nishuonimane/article/details/127247294

版权

传统奇异值分解

奇异值分解技术（Singular Value Decomposition，SVD）是一种矩阵分解方法，实际上是计算矩阵的特征值。若A是实对称矩阵，并且矩阵的A的阶次为n，则存在正交矩阵Q使得：
在这里插入图片描述其中，lamda为矩阵A的特征值。对于非实对称矩阵，式（2.24）不再适用，但存在两个正交矩阵P和V使得：
其中，T代表矩阵的转置运算，r代表矩阵U的秩，矩阵D是一个的对角矩阵，如下所示：

是矩阵U的全部奇异值，即。其中奇异值是按从大到小排列的。一般来说，信号的能量越大，它的奇异值越大。

传统奇异值分解的不足

文献1中指出传统的SVD方法可能存在两点不理想：1、尽管Householder变换本质上是平行的，但每一步的有效向量长度都会减少，从而导致效率低下；2、Sameh和Kuck的并行QR方法在数值上可能不稳定。相比之下，文献2的单边正交化方法很容易适应于并行机上的计算。所以该文介绍了一种新的奇异值的方法。

R语言代码实现

SVD_new<-function(A){
  U<-diag(nrow(A))
  V<-matrix(NA,nrow = nrow(A),ncol = ncol(A))
  squre_epm<-1*10^-30
  squre_tau<-1*10^-24
  c<-0
  squre_A<-c()
  for (i in 1:nrow(A)) {
    squre_A[i]<-norm(as.matrix(A[i,]),"F")^2
  }
  squre_delta<- squre_epm*(norm(A,"F")^2)
  for (i in 1:(nrow(A)-1)) {
    for (j in (i+1):nrow(A)) {
      if (squre_A[j]<squre_delta){
        c<-c+1
      }
      else if(squre_A[i]<squre_delta){
        C<-c()
        C<-A[i,]
        A[i,]<-A[j,]
        A[j,]<-C
        D<-c()
        D<-U[i,]
        U[i,]<-U[j,]
        U[j,]<-D
      }
      else if((sum(t(A[i,])*A[j,])/(norm(as.matrix(A[i,]),"F")*norm(as.matrix(A[j,]),"F")))^2<squre_tau){
        if(squre_A[i]<squre_A[j]){
          C<-c()
          C<-A[i,]
          A[i,]<-A[j,]
          A[j,]<-C
          D<-c()
          D<-U[i,]
          U[i,]<-U[j,]
          U[j,]<-D
        }
      }
      else{
        alpha<-2*sum(t(A[i,])*A[j,])
        beta<-squre_A[i]-squre_A[j]
        gamma<-sqrt(alpha^2+beta^2)
        if(beta>0){
          cos_fi<-sqrt((gamma+beta)/(2*gamma))
          sin_fi<-alpha/(2*gamma*cos_fi)
        }
        else{
          sin_fi<-sqrt((gamma-beta)/(2*gamma))
          cos_fi<-alpha/(2*gamma*sin_fi)
        }
        w<-cos_fi*A[i,]+sin_fi*A[j,]
        A[j,]<--sin_fi*A[i,]+cos_fi*A[j,]
        A[i,]<-w
        squre_w<-cos_fi^2*squre_A[i]+sin_fi^2*squre_A[j]+alpha*cos_fi*sin_fi
        squre_A[j]<-sin_fi^2*squre_A[i]+cos_fi^2*squre_A[j]-alpha*cos_fi*sin_fi
        squre_A[i]<-squre_w
        z<-cos_fi*U[i,]+sin_fi*U[j,]
        U[j,]<--sin_fi*U[i,]+cos_fi*U[j,]
        U[i,]<-z
      }
    }
  }
  
  
  #compute singular values
  i<-1
  sigma<-c()
  while(i<=nrow(A)&&norm(as.matrix(A[i,]),"F")>sqrt(squre_delta)){
    sigma[i]<-norm(as.matrix(A[i,]),"F")
    V[i,]<-A[i,]/norm(as.matrix(A[i,]),"F")
    i<-i+1
  }
  r<-i-1
  return(sigma)
}