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

Запустите с симметричной положительной определенной матрицы, 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' аппроксимирует A, ichol с диагональной компенсацией создает 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');