疑似逆行列を計算する Moore-Penrose アルゴリズムの Matlab 実装を探しています。
いくつかのアルゴリズムを試しましたが、これは
http://arxiv.org/ftp/arxiv/papers/0804/0804.4809.pdf
初見で良さそう。
ただし、問題は、大きな要素の場合、スケーリングが不十分な行列が生成され、一部の内部操作が失敗することです。次の手順に関係します。
L=L(:,1:r);
M=inv(L'*L);
他の SW で簡単に実装できる、より堅牢なソリューションを見つけようとしています。ご協力いただきありがとうございます。
ビルトインを使用することの何が問題になっていpinvますか？
それ以外の場合は、Octaveで使用されている実装を確認できます。Octave / MATLAB構文ではありませんが、問題なく移植できるはずです。
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;
}
これは MP 疑似逆数を計算するために書かれた [I][1] の R コードです。これは、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))
}