Использование chol факторизация матрицы симметричных коэффициентов, а затем решение линейной системы с помощью фактора Холецкого.

Создайте симметричную матрицу с положительными значениями на диагонали.

A = [1 0 1; 0 2 0; 1 0 3]

A = 3×3

     1     0     1
     0     2     0
     1     0     3

Вычислим Фактор Холецкого матрицы.

R = chol(A)

R = 3×3

    1.0000         0    1.0000
         0    1.4142         0
         0         0    1.4142

Создайте вектор для правой стороны уравнения $\mathrm{Ax}=\mathit{b}$.

b = sum(A,2);

С тех пор $\mathit{A}={\mathit{R}}^{\mathit{T}}\mathit{R}$ при разложении Холецкого линейное уравнение становится ${\mathit{R}}^{\mathit{T}}\mathit{R}\mathit{x}=\mathit{b}$. Решение для x использование оператора обратной косой черты.

x = R\(R'\b)

x = 3×1

    1.0000
    1.0000
    1.0000

Вычислите верхние и более низкие факторизации Холесского матрицы и проверьте результаты.

Создайте симметричную положительную определенную матрицу теста 6 на 6 с помощью gallery функция.

A = gallery('lehmer',6);

Вычислите фактор Холецкого, используя верхний треугольник A.

R = chol(A)

R = 6×6

    1.0000    0.5000    0.3333    0.2500    0.2000    0.1667
         0    0.8660    0.5774    0.4330    0.3464    0.2887
         0         0    0.7454    0.5590    0.4472    0.3727
         0         0         0    0.6614    0.5292    0.4410
         0         0         0         0    0.6000    0.5000
         0         0         0         0         0    0.5528

Проверьте, что верхний треугольный множитель удовлетворяет R'*R - A = 0, в пределах ошибки округления.

norm(R'*R - A)

ans = 2.5801e-16

Теперь задайте 'lower' опция вычисления фактора Холецкого с помощью нижнего треугольника A.

L = chol(A,'lower')

L = 6×6

    1.0000         0         0         0         0         0
    0.5000    0.8660         0         0         0         0
    0.3333    0.5774    0.7454         0         0         0
    0.2500    0.4330    0.5590    0.6614         0         0
    0.2000    0.3464    0.4472    0.5292    0.6000         0
    0.1667    0.2887    0.3727    0.4410    0.5000    0.5528

Проверьте, что нижний треугольный множитель удовлетворяет L*L' - A = 0, в пределах ошибки округления.

norm(L*L' - A)

ans = 2.5801e-16

Использование chol с двумя выходами, чтобы подавить ошибки, когда входная матрица не симметрично положительно определена.

Создайте матрицу 5 на 5 биномиальных коэффициентов. Эта матрица симметрично положительно определена, поэтому вычитайте 1 из последнего элемента, чтобы убедиться, что он больше не положительно определен.

A = pascal(5);
A(end) = A(end) - 1

A = 5×5

     1     1     1     1     1
     1     2     3     4     5
     1     3     6    10    15
     1     4    10    20    35
     1     5    15    35    69

Вычислим Фактор Холецкого для A. Задайте два выхода, чтобы избежать генерации ошибки, если A не симметрично положительно определено.

[R,flag] = chol(A)

R = 4×4

     1     1     1     1
     0     1     2     3
     0     0     1     3
     0     0     0     1

flag = 5

Начиная с flag является ненулевым, он задает индекс поворота, где факторизация прекращает работать. chol умеет вычислять q = flag-1 = 4 строки и столбцы правильно перед сбоем, когда они сталкиваются с измененной частью матрицы.

Проверьте, что R'*R возвращает четыре строки и столбца, которые согласуются с A(1:q,1:q).

q = flag-1;
R'*R

ans = 4×4

     1     1     1     1
     1     2     3     4
     1     3     6    10
     1     4    10    20

A(1:q,1:q)

ans = 4×4

     1     1     1     1
     1     2     3     4
     1     3     6    10
     1     4    10    20

Вычислите коэффициент Холецкого разреженной матрицы и используйте выход сочетания для создания фактора Холецкого с меньшим количеством ненулей.

Создайте разреженную положительную определенную матрицу на основе west0479 матрица.

load west0479
A = west0479;
S = A'*A;

Вычислим Фактор Холецкого матрицы двумя различными способами. Сначала задайте два выхода, а затем три выхода, чтобы разрешить изменение порядка строк и столбцов.

[R,flag] = chol(S);
[RP,flagP,P] = chol(S);

Для каждого расчета проверяйте, что flag = 0 для подтверждения успешности вычисления.

if ~flag && ~flagP
    disp('Factorizations successful.')
else
    disp('Factorizations failed.')
end

Factorizations successful.

Сравните количество ненулевых в chol(S) по сравнению с переупорядоченной матрицей chol(P'*S*P). Лучшая практика - использовать синтаксис трёх выходов chol с разреженными матрицами, поскольку переупорядочивание строк и столбцов может значительно уменьшить количество ненулей в факторе Холецкого.

subplot(1,2,1)
spy(R)
title('Nonzeros in chol(S)')
subplot(1,2,2)
spy(RP)
title('Nonzeros in chol(P''*S*P)')

Используйте 'vector' опция chol чтобы вернуть информацию о сочетании как вектор, а не как матрицу.

Создайте разреженную матрицу конечных элементов.

S = gallery('wathen',10,10);
spy(S)

Вычислите фактор Холецкого для матрицы и задайте 'vector' опция для возврата вектора сочетания p.

[R,flag,p] = chol(S,'vector');

Проверьте, что flag = 0, что указывает на успешность вычисления.

if ~flag
    disp('Factorization successful.')
else
    disp('Factorization failed.')
end

Factorization successful.

Проверьте, что S(p,p) = R'*R, в пределах ошибки округления.

norm(S(p,p) - R'*R,'fro')

ans = 2.1039e-13