Matlab: реализация псевдообратного алгоритма Мура-Пенроуза

Я ищу реализацию Matlab алгоритма Мура-Пенроуза, вычисляющего псевдообратную матрицу.

Я пробовал несколько алгоритмов, этот

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

показался хорошим на первый взгляд.

Однако проблема в том, что для больших элементов он создает плохо масштабируемые матрицы, и некоторые внутренние операции не выполняются. Это касается следующих шагов:

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

Я пытаюсь найти более надежное решение, которое легко реализовать в моем другом ПО. Спасибо за вашу помощь.


matlab matrix-inverse
person justik    schedule 11.11.2012    source источник
comment
Вы можете задать этот вопрос либо на math.stackexchange.com, либо на dsp.stackexchange.com   -  person Ali    schedule 11.11.2012

Ответы (2)


arrow_upward
4
arrow_downward

Я повторно реализовал один на C#, используя матричную библиотеку Mapack от Lutz Roeder. Возможно, эта или 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;
}
person Frank Hileman    schedule 08.12.2012

arrow_upward
4
arrow_downward

Что не так с использованием встроенного pinv?

В противном случае вы можете взглянуть на реализацию, используемую в Octave. Это не в синтаксисе Octave/MATLAB, но я думаю, вы сможете без особых проблем перенести его.

person Egon    schedule 11.11.2012
comment
Да, с пинвом все в порядке. Но я хотел бы использовать код в другом ПО, написанном на другом языке. - person justik; 11.11.2012
comment
Я думаю, что псевдоинверсия должна быть доступна практически для любого приличного языка программирования (например, с использованием библиотеки LAPACK). В общем, я бы не рекомендовал самостоятельно реализовывать численные алгоритмы для всего, что должно быть надежным (если, конечно, вы не знаете, что делаете). - person Egon; 12.11.2012