标签: svd

我正在尝试解决一组形式为Ax = 0的方程式.A是已知的6x6矩阵,我使用SVD编写了下面的代码,以获得在某种程度上起作用的向量x.答案大致正确,但不足以对我有用,我怎样才能提高计算的精确度？降低低于1.e-4的eps会导致功能失败.

from numpy.linalg import *
from numpy import *

A = matrix([[0.624010149127497 ,0.020915658603923 ,0.838082638087629 ,62.0778180312547 ,-0.336 ,0],
[0.669649399820597 ,0.344105317421833 ,0.0543868015800246 ,49.0194290212841 ,-0.267 ,0],
[0.473153758252885 ,0.366893577716959 ,0.924972565581684 ,186.071352614705 ,-1 ,0],
[0.0759305208803158 ,0.356365401030535 ,0.126682113674883 ,175.292109352674 ,0 ,-5.201],
[0.91160934274653 ,0.32447818779582 ,0.741382053883291 ,0.11536775372698 ,0 ,-0.034],
[0.480860406786873 ,0.903499596111067 ,0.542581424762866 ,32.782593418975 ,0 ,-1]])

def null(A, eps=1e-3):
  u,s,vh = svd(A,full_matrices=1,compute_uv=1)
  null_space = compress(s <= eps, vh, axis=0)
  return null_space.T

NS = null(A)
print "Null space equals ",NS,"\n"
print dot(A,NS)

与sklearn的PCA不同,TruncatedSVD的解释方差比率不是降序.我查看了源代码,似乎他们使用不同的方式计算解释的方差比:

TruncatedSVD:

U, Sigma, VT = randomized_svd(X, self.n_components,
                              n_iter=self.n_iter,
                              random_state=random_state)
X_transformed = np.dot(U, np.diag(Sigma))
self.explained_variance_ = exp_var = np.var(X_transformed, axis=0)
if sp.issparse(X):
    _, full_var = mean_variance_axis(X, axis=0)
    full_var = full_var.sum()
else:
    full_var = np.var(X, axis=0).sum()
self.explained_variance_ratio_ = exp_var / full_var

PCA:

U, S, V = linalg.svd(X, full_matrices=False)
explained_variance_ = (S ** 2) / n_samples
explained_variance_ratio_ = (explained_variance_ /
                             explained_variance_.sum())

PCA使用sigma直接计算explain_variance,由于sigma按降序排列,explain_variance也按降序排列.另一方面,TruncatedSVD使用变换矩阵的列的方差来计算explain_variance,因此方差不一定按降序排列.

这是否意味着我需要explained_variance_ratio从头TruncatedSVD开始排序才能找到前k个主要组件？

我的MatrixR 稀疏,显然对我来说太大as.matrix()了(虽然它也不是超大).有as.matrix()问题的调用是在svd()函数内部,所以我想知道是否有人知道不需要先转换为密集矩阵的SVD的不同实现.

我试图在R中计算P任意N x J矩阵的投影矩阵S:

P = S (S'S) ^ -1 S'

我一直试图用以下功能执行此操作:

P <- function(S){
  output <- S %*% solve(t(S) %*% S) %*% t(S)
  return(output)
}

但是当我使用它时,我得到的错误看起来像这样:

# Error in solve.default(t(S) %*% S, t(S), tol = 1e-07) : 
#  system is computationally singular: reciprocal condition number = 2.26005e-28

我认为这是数值下溢和/或不稳定的结果在许多地方一样讨论R-帮助和这里,但我没有足够的经验使用SVD或QR分解来解决这个问题,要不然就把这个现有的代码到行动.我也尝试了建议的代码,即将solve写成一个系统:

output <- S %*% solve (t(S) %*% S, t(S), tol=1e-7)

但它仍然无效.任何建议,将不胜感激.

我很确定我的矩阵应该是可逆的并且没有任何共线性,只是因为我尝试用正交虚拟变量矩阵进行测试,但它仍然不起作用.

另外,我想将它应用于相当大的矩阵,所以我正在寻找一个简洁的通用解决方案.

尝试在Python中计算SVD以找到频谱中最重要的元素,并创建一个仅包含最重要部分的矩阵.

在python我有:

u,s,v = linalg.svd(Pxx, full_matrices=True)

这样可以得到3个矩阵; 其中"s"包含与u,v相对应的大小.

为了构造一个包含信号的所有重要部分的新矩阵,我需要捕获"s"中的最高值并将它们与"u"和"v"中的列匹配,得到的矩阵应该给我数据中最重要的部分.

问题是我不知道如何在Python中执行此操作,例如,如何在"s"中找到最高数字并选择"u"和"v"中的列以创建新矩阵？

(我是Python和numpy的新手)所以任何帮助都将非常感激

编辑:

import wave, struct, numpy as np, matplotlib.mlab as mlab, pylab as pl
from scipy import linalg, mat, dot;
def wavToArr(wavefile):
    w = wave.open(wavefile,"rb")
    p = w.getparams()
    s = w.readframes(p[3])
    w.close()
    sd = np.fromstring(s, np.int16)
    return sd,p

def wavToSpec(wavefile,log=False,norm=False):
    wavArr,wavParams = wavToArr(wavefile)
    print wavParams
    return  mlab.specgram(wavArr, NFFT=256,Fs=wavParams[2],detrend=mlab.detrend_mean,window=mlab.window_hanning,noverlap=128,sides='onesided',scale_by_freq=True)

wavArr,wavParams = wavToArr("wavBat1.wav")

Pxx, freqs, bins = wavToSpec("wavBat1.wav")
Pxx += 0.0001

U, s, Vh = linalg.svd(Pxx, full_matrices=True)
assert np.allclose(Pxx, np.dot(U, np.dot(np.diag(s), Vh))) …

我想在PHP中实现奇异值分解(SVD).我知道有几个外部库可以为我做这个.但是我有两个关于PHP的问题:1)你认为用PHP编写SVD是否可能和/或合理？2)如果(1)是肯定的:你能帮我用PHP编写代码吗？

我已经自己编写了SVD的一些部分.这是我对行动过程做出评论的代码.此代码的某些部分并不完全正确.

如果你能帮助我会很棒.非常感谢你提前!

我正在使用Octave和R来使用简单的矩阵计算SVD并得到两个不同的答案!代码如下:

[R

> a<-matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,1,1,1), 9, 4)
> a
 [,1] [,2] [,3] [,4]
 [1,]    1    1    0    0
 [2,]    1    1    0    0
 [3,]    1    1    0    0
 [4,]    1    0    1    0
 [5,]    1    0    1    0
 [6,]    1    0    1    0
 [7,]    1    0    0    1
 [8,]    1    0    0    1
 [9,]    1    0    0    1
> a.svd <- svd(a)
> a.svd$d
[1] 3.464102e+00 1.732051e+00 1.732051e+00 1.922963e-16
> a.svd$u
       [,1]       [,2]          [,3]          [,4]
 [1,] -0.3333333  0.4714045 -1.741269e-16  7.760882e-01
 [2,] …

更新:我已经修改了Optimize,Eigen和Solve方法来反映变化.所有现在都返回"相同"向量,允许机器精度. 我仍然对Eigen方法感到困惑.具体如何/为什么我选择特征向量切片没有意义.这只是试验和错误,直到正常情况与其他解决方案相匹配.如果有人能够纠正/解释我真正应该做的事情,或者为什么我做了什么工作我会很感激..

感谢 亚历山大克莱默解释为什么我采取切片,只能选择一个正确的答案

我有一个深度图像.我想计算深度图像中像素的粗糙表面法线.我考虑周围像素,在最简单的情况下是3x3矩阵,并将平面拟合到这些点,并计算该平面的法向单位矢量.

听起来很简单,但最好先考虑平面拟合算法.搜索SO和其他各种网站,我看到使用最小二乘法,单值分解,特征向量/值等的方法.

虽然我不完全理解数学,但我已经能够使各种片段/示例起作用.我遇到的问题是,我对每种方法都有不同的答案.我期待各种答案是相似的(不完全),但它们似乎有很大的不同.也许有些方法不适合我的数据,但不确定为什么我会得到不同的结果.有什么想法吗？

这是代码的更新输出:

LTSQ:   [ -8.10792259e-17   7.07106781e-01  -7.07106781e-01]
SVD:    [ 0.                0.70710678      -0.70710678]
Eigen:  [ 0.                0.70710678      -0.70710678]
Solve:  [ 0.                0.70710678       0.70710678]
Optim:  [ -1.56069661e-09   7.07106781e-01   7.07106782e-01]

以下代码实现了五种不同的方法来计算平面的曲面法线.算法/代码来自互联网上的各种论坛.

import numpy as np
import scipy.optimize

def fitPLaneLTSQ(XYZ):
    # Fits a plane to a point cloud, 
    # Where Z = aX + bY + c        ----Eqn #1
    # Rearanging Eqn1: aX + bY -Z +c =0
    # Gives normal (a,b,-1)
    # …

SciPy和Numpy都内置了奇异值分解(SVD)函数.命令基本上scipy.linalg.svd和numpy.linalg.svd.这两者有什么区别？他们中的任何一个比另一个好吗？

我校准了相机并找到了内在参数(K).我也计算了基本矩阵(F).

现在E = K_T*F*K. 到现在为止还挺好.

现在我们将基本矩阵(E)传递给SVD,使用分解值(U,W,V)来提取旋转和平移:

 essentialMatrix = K.Transpose().Mul(fund).Mul(K);
 CvInvoke.cvSVD(essentialMatrix, wMatrix, uMatrix, vMatrix, Emgu.CV.CvEnum.SVD_TYPE.CV_SVD_DEFAULT);

**问题)此时,已经提出了两种方法,并且让我感到困惑的是哪一种方法确实给出了正确的答案 - 特别是对于翻译:

首先输入链接描述,作者建议计算R,T如下:

但在第二种方法[ http://isit.u-clermont1.fr/~ab/Classes/DIKU-3DCV2/Handouts/Lecture16.pdf]中,作者为T提供了另一个公式,即+ U,-U,如下所示:

我正在使用openCv库在C#.Net上实现它.谁知道哪个翻译公式是正确的？