В этом примере создается неполная факторизация Холеского.

Начните с симметричной положительной определенной матрицы, 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));

Создание неполной факторизации Cholesky в качестве предварительного условия для 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). Поскольку эта матрица является дискретизированной лапласианской, однако, использование модифицированной неполной Cholesky может создать лучшее предварительное условие. Модифицированная неполная факторизация Холеского конструирует приближённую факторизацию, сохраняющую действие оператора на вектор константы. То есть 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)');

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

Также можно попробовать неполные факторизации Cholesky с понижением порога. Следующий график показывает сходимость 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');

Обратите внимание на неполное предварительное кондиционирующее устройство Cholesky, построенное с допуском падения 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');