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

За 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 обновил секция алгоритма на странице документа mldivide с некоторыми хорошими блок-схемами. Увидеть Вот а также Вот (полные и редкие случаи).

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

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

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

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

x = pinv(A)*b

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

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

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

Более того, оператор обратной косой черты также работает на gpuArrayв этом случае он опирается на cuBLAS а также МАГМА выполнить на GPU.

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

Это довольно сложный заказ, если вы хотите реализовать все это сами 🙂