2012-11-11 81 views
3

我正在尋找一個計算僞逆矩陣的Moore-Penrose算法的matlab實現。Matlab:Moore-Penrose僞逆算法的實現

我試了幾個algoithm,這一個

http://arxiv.org/ftp/arxiv/papers/0804/0804.4809.pdf

一見鍾情出現良好。

但是,它的問題是,對於大型元素,它會產生嚴重縮放的矩陣和一些內部操作失敗。它涉及以下步驟:

L=L(:,1:r); 
M=inv(L'*L); 

我試圖找到一個更強大的解決方案,它是我在其他SW容易實現。謝謝你的幫助。

+1

你可能想利用這個問題要麼http://math.stackexchange.com/或http://dsp.stackexchange.com/ – Ali

回答

4

我使用Lutz Roeder的Mapack矩陣庫在C#中重新實現了一個。也許這個,或者Java版本,對你有用。

/// <summary> 
/// The difference between 1 and the smallest exactly representable number 
/// greater than one. Gives an upper bound on the relative error due to 
/// rounding of floating point numbers. 
/// </summary> 
const double MACHEPS = 2E-16; 

// NOTE: Code for pseudoinverse is from: 
// http://the-lost-beauty.blogspot.com/2009/04/moore-penrose-pseudoinverse-in-jama.html 

/// <summary> 
/// Computes the Moore–Penrose pseudoinverse using the SVD method. 
/// Modified version of the original implementation by Kim van der Linde. 
/// </summary> 
/// <param name="x"></param> 
/// <returns>The pseudoinverse.</returns> 
public static Matrix MoorePenrosePsuedoinverse(Matrix x) 
{ 
    if (x.Columns > x.Rows) 
     return MoorePenrosePsuedoinverse(x.Transpose()).Transpose(); 
    SingularValueDecomposition svdX = new SingularValueDecomposition(x); 
    if (svdX.Rank < 1) 
     return null; 
    double[] singularValues = svdX.Diagonal; 
    double tol = Math.Max(x.Columns, x.Rows) * singularValues[0] * MACHEPS; 
    double[] singularValueReciprocals = new double[singularValues.Length]; 
    for (int i = 0; i < singularValues.Length; ++i) 
     singularValueReciprocals[i] = Math.Abs(singularValues[i]) < tol ? 0 : (1.0/singularValues[i]); 
    Matrix u = svdX.GetU(); 
    Matrix v = svdX.GetV(); 
    int min = Math.Min(x.Columns, u.Columns); 
    Matrix inverse = new Matrix(x.Columns, x.Rows); 
    for (int i = 0; i < x.Columns; i++) 
     for (int j = 0; j < u.Rows; j++) 
      for (int k = 0; k < min; k++) 
       inverse[i, j] += v[i, k] * singularValueReciprocals[k] * u[j, k]; 
    return inverse; 
} 
4

使用內置pinv有什麼問題?您可以看看implementation used in Octave。它不在Octave/MATLAB語法中,但我想你應該能夠移植它而沒有太多問題。

+0

是,PINV是OK 。但我想在另一個用不同語言編寫的軟件中使用代碼。 – justik

+0

我認爲僞逆應該適用於幾乎任何體面的編程語言(例如使用LAPACK庫)。一般來說,我不建議你自己實施數值算法,以確保任何應該可靠的數據(除非你知道你在做什麼)。 – Egon

1

這裏是R代碼[I] [1]已經寫出來計算M-P pseudoinverse。我認爲這很簡單,可以轉換成matlab代碼。

pinv<-function(H){ 
 
    x=t(H) %*% H 
 
    s=svd(x) 
 
    xp=s$d 
 
    for (i in 1:length(xp)){ 
 
    if (xp[i] != 0){ 
 
     xp[i]=1/xp[i] 
 
    } 
 
    else{ 
 
     xp[i]=0 
 
    } 
 
    } 
 
    return(s$u %*% diag(xp) %*% t(s$v) %*% t(H)) 
 
}
[1]: http://hamedhaseli.webs.com/downloads