我正在寻找一个Matlab实现的Moore-Penrose算法,用于计算伪逆矩阵。我尝试了几种算法,这个链接看起来很不错。然而,问题在于,对于大元素,它会产生缩放不良的矩阵,一些内部操作会失败。这涉及以下步骤:
我正在尝试寻找一种更加健壮的解决方案,该方案可以轻松地在我的其他软件中实现。感谢您的帮助。
L=L(:,1:r);
M=inv(L'*L);
我正在尝试寻找一种更加健壮的解决方案,该方案可以轻松地在我的其他软件中实现。感谢您的帮助。
L=L(:,1:r);
M=inv(L'*L);
我使用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;
}
使用内置的pinv
有什么问题吗?
否则,您可以查看Octave中使用的实现方式。它不是用Octave / MATLAB语法编写的,但我想您应该能够轻松地进行移植。
Here is the R code [I][1] have written to compute M-P pseudoinverse. I think that is simple enough to be translated into matlab code.
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))
}