Этот пример генерирует неполную факторизацию Холесского.

Начните с симметричной положительной определенной матрицы, A:

N = 100;
A = delsq(numgrid('S',N));

A двумерный, дискретный отрицательный Лапласиан с пятью точками на 100 100 квадратной сетке с граничными условиями Дирихле. Размер A 98*98 = 9604 (не 10000, когда границы сетки используются, чтобы наложить условия Дирихле).

Неполная факторизация Холесского без заливок является факторизацией, которая содержит только ненули в той же позиции A содержит ненули. Эта факторизация является чрезвычайно дешевой, чтобы вычислить. Несмотря на то, что продукт L*L' обычно очень отличается от A, продукт L*L' будет совпадать с A на его шаблоне до округления.

L = ichol(A);
norm(A-L*L','fro')./norm(A,'fro')

ans = 0.0916

norm(A-(L*L').*spones(A),'fro')./norm(A,'fro')

ans = 4.9606e-17

ichol май также использоваться, чтобы сгенерировать неполные факторизации Холесского с пороговым отбрасыванием. Когда допуск отбрасывания уменьшается, фактор имеет тенденцию становиться более плотным и продукт L*L' имеет тенденцию быть лучшим приближением A. Следующие графики показывают относительную погрешность неполной факторизации, построенной против допуска отбрасывания, а также отношения плотности неполных факторов к плотности полного Фактора Холесского.

n = size(A,1);
ntols = 20;
droptol = logspace(-8,0,ntols);
nrm = zeros(1,ntols);
nz = zeros(1,ntols);
nzComplete = nnz(chol(A,'lower'));
for k = 1:ntols
    L = ichol(A,struct('type','ict','droptol',droptol(k)));
    nz(k) = nnz(L);
    nrm(k) = norm(A-L*L','fro')./norm(A,'fro');
end
figure
loglog(droptol,nrm,'LineWidth',2)
title('Drop tolerance vs norm(A-L*L'',''fro'')./norm(A,''fro'')')

figure
semilogx(droptol,nz./nzComplete,'LineWidth',2)
title('Drop tolerance vs fill ratio ichol/chol')

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

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

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

N = 100;
A = delsq(numgrid('S',N));

Создайте неполную факторизацию Холесского как предварительный формирователь для pcg. Используйте постоянный вектор в качестве правой стороны. Как базовая линия, выполните pcg без предварительного формирователя.

b = ones(size(A,1),1);
tol = 1e-6;
maxit = 100;
[x0,fl0,rr0,it0,rv0] = pcg(A,b,tol,maxit);

Обратите внимание на то, что fl0 = 1 указание, что pcg не управлял относительной невязкой к требуемому допуску в максимальных позволенных итерациях. Попробуйте неполную факторизацию Холесского без заливок как предварительный формирователь.

L1 = ichol(A);
[x1,fl1,rr1,it1,rv1] = pcg(A,b,tol,maxit,L1,L1');

fl1 = 0, указание, что pcg сходившийся к требуемому допуску и сделал так в 59 итерациях (значение it1). Поскольку этой матрицей является дискретизированный Лапласиан, однако, использование изменило неполный Холесский, может создать лучший предварительный формирователь. Модифицированная неполная факторизация Холесского создает аппроксимированную факторизацию, которая сохраняет действие оператора на постоянном векторе. Таким образом, norm(A*e-L*(L'*e)) будет приблизительно нуль для e = ones(size(A,2),1) даже при том, что norm(A-L*L','fro')/norm(A,'fro') не близко к нулю. Не необходимо задать тип для этого синтаксиса начиная с nofill значение по умолчанию, но это - хорошая практика.

opts.type = 'nofill';
opts.michol = 'on';
L2 = ichol(A,opts);
e = ones(size(A,2),1);
norm(A*e-L2*(L2'*e))

ans = 3.7983e-14

[x2,fl2,rr2,it2,rv2] = pcg(A,b,tol,maxit,L2,L2');

pcg сходится (fl2 = 0) но только в 38 итерациях. Графический вывод всех трех историй сходимости показывает сходимость.

semilogy(0:maxit,rv0./norm(b),'b.');
hold on
semilogy(0:it1,rv1./norm(b),'r.');
semilogy(0:it2,rv2./norm(b),'k.');
legend('No Preconditioner','IC(0)','MIC(0)');

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

Можно также попробовать неполные факторизации Холесского пороговым отбрасыванием. Следующий график показывает сходимость pcg с предварительными формирователями, созданными с различными допусками отбрасывания.

L3 = ichol(A, struct('type','ict','droptol',1e-1));
[x3,fl3,rr3,it3,rv3] = pcg(A,b,tol,maxit,L3,L3');
L4 = ichol(A, struct('type','ict','droptol',1e-2));
[x4,fl4,rr4,it4,rv4] = pcg(A,b,tol,maxit,L4,L4');
L5 = ichol(A, struct('type','ict','droptol',1e-3));
[x5,fl5,rr5,it5,rv5] = pcg(A,b,tol,maxit,L5,L5'); 

figure
semilogy(0:maxit,rv0./norm(b),'b-','linewidth',2);
hold on
semilogy(0:it3,rv3./norm(b),'b-.','linewidth',2);
semilogy(0:it4,rv4./norm(b),'b--','linewidth',2);
semilogy(0:it5,rv5./norm(b),'b:','linewidth',2);
legend('No Preconditioner','ICT(1e-1)','ICT(1e-2)', ...
   'ICT(1e-3)','Location','SouthEast');

Отметьте неполный предварительный формирователь Холесского, созданный с допуском отбрасывания 1e-2 обозначается как ICT(1e-2).

Как с заполнением нулями неполного Холесского, пороговая факторизация отбрасывания может извлечь выгоду из модификации (т.е. opts.michol = 'on') поскольку матрица является результатом эллиптического дифференциального уравнения с частными производными. Как с MIC(0), модифицированный пороговый понижающийся неполный Холесский сохранит действие предварительного формирователя на постоянных векторах, который является norm(A*e-L*(L'*e)) будет приблизительно нуль.

Этот пример иллюстрирует использование diagcomp опция ichol.

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

S = rng('default');
A = sprandsym(1000,1e-2,1e-4,1);
rng(S);
b = full(sum(A,2));
[x0,fl0,rr0,it0,rv0] = pcg(A,b,1e-6,100);

Поскольку сходимость не достигнута, попытайтесь создать неполный предварительный формирователь Холесского.

L = ichol(A,struct('type','ict','droptol',1e-3));

Error using ichol
Encountered nonpositive pivot.

Если ichol ломается как выше, можно использовать diagcomp опция, чтобы создать переключенную неполную факторизацию Холесского. Таким образом, вместо того, чтобы создать L таким образом, что L*L' аппроксимирует Aichol с диагональной компенсацией создает L таким образом, что L*L' аппроксимирует M = A + alpha*diag(diag(A)) явным образом не формируя M. Когда неполные факторизации всегда существуют для по диагонали доминирующих матриц, alpha как могут находить, делает M по диагонали доминирующий.

alpha = max(sum(abs(A),2)./diag(A))-2

alpha = 62.9341

L1 = ichol(A, struct('type','ict','droptol',1e-3,'diagcomp',alpha));
[x1,fl1,rr1,it1,rv1] = pcg(A,b,1e-6,100,L1,L1');

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

alpha = .1;
L2 = ichol(A, struct('type','ict','droptol',1e-3,'diagcomp',alpha));
[x2,fl2,rr2,it2,rv2] = pcg(A,b,1e-6,100,L2,L2');

Теперь pcg сходится и график может показать истории сходимости с каждым предварительным формирователем.

semilogy(0:100,rv0./norm(b),'b.');
hold on;
semilogy(0:100,rv1./norm(b),'r.');
semilogy(0:it2,rv2./norm(b),'k.');
legend('No Preconditioner','\alpha \approx 63','\alpha = .1');
xlabel('Iteration Number');
ylabel('Relative Residual');