ichol

Неполная факторизация Холесского

Синтаксис

L = ichol(A)
L = ichol(A,opts)

Описание

L = ichol(A) выполняет неполную факторизацию Холесского A с, заполняют нулями.

L = ichol(A,opts) выполняет неполную факторизацию Холесского A с опциями, заданными opts.

По умолчанию ichol ссылается на более низкий треугольник A и производит нижние треугольные факторы.

Входные параметры

A

Разреженная матрица

opts

Структура максимум с пятью полями:

Имя поляСводные данныеОписание
typeТип факторизации

Указывает который разновидность неполного Холесского выполнить. Допустимыми значениями этого поля является 'nofill' и 'ict'. Вариант 'nofill' выполняет incomplete Cholesky with zero-fill (IC(0)). Вариант 'ict' performs incomplete Cholesky with threshold dropping (ICT). Значением по умолчанию является 'nofill'.

droptol Пропустите допуск, когда типом будет 'ict'

Неотрицательный скаляр, используемый в качестве допуска отбрасывания при выполнении ICT. Элементы, которые меньше в значении, чем локальный допуск отбрасывания, исключены из получившегося фактора за исключением диагонального элемента, который никогда не уронен. Локальным допуском отбрасывания на шаге j факторизации является norm(A(j:end,j),1)*droptol. 'droptol' проигнорирован, если 'type' является 'nofill'. Значением по умолчанию является 0.

micholУказывает, выполнить ли измененный неполный Холесский

Указывает, выполняется ли modified incomplete Cholesky (MIC). Полем может быть 'on' или 'off'. При выполнении MIC диагональ компенсируется уроненные элементы, чтобы осуществить отношение A*e = L*L'*e где   e = ones(size(A,2),1). Значением по умолчанию является 'off'.

diagcompВыполните компенсированный неполный Холесский с заданным коэффициентом

Действительный неотрицательный скаляр, используемый в качестве глобальной диагонали, переключает alpha в формировании неполного Фактора Холесского. Таким образом, вместо того, чтобы выполнить неполный Холесский на A, формируется факторизация A + alpha*diag(diag(A)). Значением по умолчанию является 0.

shapeОпределяет, на какой треугольник ссылаются и возвращают

Допустимыми значениями является 'upper' и 'lower'. Если 'upper' задан, только на верхний треугольник A ссылаются, и R создается таким образом, что A аппроксимирован R'*R. Если 'lower' задан, только на более низкий треугольник A ссылаются, и L создается таким образом, что A аппроксимирован L*L'. Значением по умолчанию является 'lower'.

Примеры

свернуть все

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

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

Советы

  • Фактор, данный этой стандартной программой, может быть полезен как предварительный формирователь для системы линейных уравнений, решаемых итерационными методами, такими как pcg или minres.

  • ichol работает только на разреженные квадратные матрицы

Ссылки

[1] Саад, Yousef. “Предварительное создание условий методов”. Итерационные методы для разреженных линейных систем. PWS Publishing Company, 1996.

[2] Мантойфель, T.A. “Неполный метод факторизации для положительных определенных линейных систем”. Математика. Comput. 34, 473–497, 1980.

Смотрите также

| | |

Для просмотра документации необходимо авторизоваться на сайте