Используйте 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.1194e-13