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

Начнем с симметричной положительно определенной матрицы, 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 по умолчанию, но это хорошая практика.

options.type = 'nofill';
options.michol = 'on';
L2 = ichol(A,options);
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).

Как и в случае неполного факторизации Холецкого с нулевым заполнением, факторизация порогового падения может быть полезной от модификации (то есть options.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');