Использовать 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{=}}^{\mathit{}}\mathit{}$RTR с разложением Холеского, линейное уравнение ${\mathit{становится}}^{\mathit{}}\mathit{}\mathit{RTR}\mathit{x}$= 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