Как реализовать Matlab mldivide (a.k.a. оператор обратной косой черты "\" )

В настоящее время я пытаюсь разработать небольшую матрично-ориентированную математическую библиотеку (я использую Eigen 3 для матричных структур данных и операций) и я хотел бы реализовать некоторые удобные функции Matlab, такие как широко используемый оператор обратной косой черты (который эквивалентен mldivide), чтобы вычислить решение линейных систем (выраженное в матричной форме).

Есть ли хорошее подробное объяснение того, как это можно достичь? (Я уже реализовал псевдо-обратную функцию Moore-Penrose pinv с классической декомпозицией SVD, но я где-то читал, что A\b не всегда pinv(A)*b, по крайней мере, Matalb не просто делает это)

Спасибо

Ответ 1

Для x = A\b оператор обратная косая черта включает в себя несколько алгоритмов для обработки различных типов ввода матрицы. Таким образом, матрица A диагностируется и выбирается путь выполнения в соответствии с ее характеристиками.

Следующая страница описывает в псевдокоде, когда A является полной матрицей:

if size(A,1) == size(A,2)         % A is square
    if isequal(A,tril(A))         % A is lower triangular
        x = A \ b;                % This is a simple forward substitution on b
    elseif isequal(A,triu(A))     % A is upper triangular
        x = A \ b;                % This is a simple backward substitution on b
    else
        if isequal(A,A')          % A is symmetric
            [R,p] = chol(A);
            if (p == 0)           % A is symmetric positive definite
                x = R \ (R' \ b); % a forward and a backward substitution
                return
            end
        end
        [L,U,P] = lu(A);          % general, square A
        x = U \ (L \ (P*b));      % a forward and a backward substitution
    end
else                              % A is rectangular
    [Q,R] = qr(A);
    x = R \ (Q' * b);
end

Для неквадратных матриц используется QR-декомпозиция. Для квадратных треугольных матриц он выполняет простую форвардную/обратную замену. Для квадратных симметричных положительно определенных матриц используется декомпозиция Холецкого. В противном случае декомпозиция LU > используется для общих квадратных матриц.

Обновление: MathWorks обновила раздел алгоритм на странице doc mldivide с некоторыми хорошими блок-схемами. См. здесь и здесь (полные и редкие случаи).

Все эти алгоритмы имеют соответствующие методы в LAPACK, и на самом деле это, вероятно, то, что делает MATLAB (обратите внимание, что последние версии корабля MATLAB с оптимизированным Intel MKL).

Причина наличия разных методов заключается в том, что он пытается использовать наиболее специфический алгоритм для решения системы уравнений, которая использует все характеристики матрицы коэффициентов (либо потому, что она будет более быстрой, либо более численной стабильностью). Поэтому вы, безусловно, можете использовать общий решатель, но он не будет наиболее эффективным.

Фактически, если вы знаете, что такое A как раньше, вы можете пропустить дополнительный процесс тестирования, вызвав linsolve и указав параметры напрямую.

Если A является прямоугольным или сингулярным, вы также можете использовать PINV, чтобы найти минимальное решение наименьших квадратов (реализовано с помощью Разложение SVD):

x = pinv(A)*b

Все вышесказанное относится к плотным матрицам, редкие матрицы - совсем другая история. Обычно в таких случаях используются итерационные решатели. Я считаю, что MATLAB использует UMFPACK и другие связанные библиотеки из пакета SuiteSpase для прямых решателей.

При работе с разреженными матрицами вы можете включить диагностическую информацию и просмотреть выполненные тесты и алгоритмы, используя spparms:

spparms('spumoni',2)
x = A\b;

Что еще, оператор обратной косой черты также работает на gpuArray, и в этом случае он полагается на cuBLAS и MAGMA для выполнения на графическом процессоре.

Он также реализован для распределенных массивов, который работает в распределенной вычислительной среде (работа разделена между кластерами компьютеров, где каждый рабочий имеет только часть массива, возможно где вся матрица не может быть сохранена в памяти сразу). В основе реализации лежит ScaLAPACK.

Это довольно высокий порядок, если вы хотите реализовать все это самостоятельно:)