symmlq

Решите систему линейных уравнений — симметричный метод LQ

Описание

пример

x = symmlq(A,b) попытки решить систему линейных уравнений A*x = b для x использование Симметричного Метода LQ. Когда попытка успешна, symmlq отображает сообщение, чтобы подтвердить сходимость. Если symmlq сбои, чтобы сходиться после максимального количества итераций или остановов по любой причине, это отображает диагностическое сообщение, которое включает относительный остаточный norm(b-A*x)/norm(b) и номер итерации, в который остановленный метод.

пример

x = symmlq(A,b,tol) задает допуск к методу. Допуском по умолчанию является 1e-6.

пример

x = symmlq(A,b,tol,maxit) задает максимальное количество итераций, чтобы использовать. symmlq отображает диагностическое сообщение, если ему не удается сходиться в maxit итерации.

пример

x = symmlq(A,b,tol,maxit,M) задает матрицу перед формирователем M и вычисляет x путем эффективного решения системы H1AHTy=H1b для y, где y=HTx и H=M1/2=(M1M2)1/2. Алгоритм не формирует H явным образом. Используя предварительный формирователь матрица может улучшить числовые свойства проблемы и КПД вычисления.

пример

x = symmlq(A,b,tol,maxit,M1,M2) задает факторы матрицы перед формирователем M таким образом, что M = M1*M2.

пример

x = symmlq(A,b,tol,maxit,M1,M2,x0) задает исходное предположение для вектора решения x. Значением по умолчанию является нулевой вектор.

пример

[x,flag] = symmlq(___) возвращает флаг, который задает, сходился ли алгоритм успешно. Когда flag = 0, сходимость была успешна. Можно использовать этот выходной синтаксис с любой из предыдущих комбинаций входных аргументов. Когда вы задаете flag вывод , symmlq не отображает диагностических сообщений.

пример

[x,flag,relres] = symmlq(___) также возвращает остаточную ошибку в вычисленном решении. Если flag 0, затем relres <= tol.

пример

[x,flag,relres,iter] = symmlq(___) также возвращает номер итерации iter в котором x был вычислен.

пример

[x,flag,relres,iter,resvec] = symmlq(___) также возвращает вектор из нормы невязки в каждой итерации, включая первый остаточный norm(b-A*x0).

пример

[x,flag,relres,iter,resvec,resveccg] = symmlq(___) также возвращает вектор из норм невязки методов сопряженных градиентов в каждой итерации.

Примеры

свернуть все

Решите квадратную линейную систему с помощью symmlq с настройками по умолчанию, и затем настраивают допуск и количество итераций, используемых в процессе решения.

Создайте разреженный трехдиагональный матричный A как матрица коэффициентов. Используйте плотные суммы строки A как вектор для правой стороны Ax=b так, чтобы решение x как ожидают, будет вектором из единиц.

n = 400; 
on = ones(n,1); 
A = spdiags([-2*on 4*on -2*on],-1:1,n,n);
b = full(sum(A,2));

Решить Ax=b использование symmlq. Выходное отображение включает значение относительной остаточной ошибки b-Axb.

x = symmlq(A,b);
symmlq stopped at iteration 20 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 20) has relative residual 0.045.

symmlq по умолчанию использование 20 итераций и допуск 1e-6, и алгоритм не может сходиться в тех 20 итерациях для этой матрицы. Поскольку невязка находится порядка 1e-2, это - хороший индикатор, что необходимо больше итераций. Также можно использовать больший допуск, чтобы облегчить для алгоритма сходиться.

Решите систему снова с помощью допуска 1e-4 и 250 итераций.

x = symmlq(A,b,1e-4,250);
symmlq converged at iteration 199 to a solution with relative residual 1.5e-14.

Исследуйте эффект использования матрицы перед формирователем с symmlq решить линейную систему.

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

A = delsq(numgrid('S',102));

Задайте b так, чтобы истинное решение Ax=b вектор из всех единиц.

b = sum(A,2);

Установите погрешность и максимальное количество итераций.

tol = 1e-12;
maxit = 100;

Используйте symmlq найти решение в требуемом допуске и количестве итераций. Задайте шесть выходных параметров, чтобы возвратить информацию о процессе решения:

  • x вычисленное решение A*x = b.

  • fl0 флаг, указывающий, сходился ли алгоритм.

  • rr0 относительная невязка вычисленного ответа x.

  • it0 номер итерации когда x был вычислен.

  • rv0 вектор из остаточной истории для b-Ax.

  • rvcg0 вектор из истории невязки метода сопряженных градиентов для ATAx-ATb.

[x,fl0,rr0,it0,rv0,rvcg0] = symmlq(A,b,tol,maxit);
fl0
fl0 = 1
rr0
rr0 = 0.0031
it0
it0 = 100

fl0 1 потому что symmlq не сходится к требуемому допуску 1e-12 в требуемых 100 итерациях.

Чтобы помочь со сходимостью, можно задать матрицу перед формирователем. Начиная с A симметрично, используйте ichol сгенерировать предварительный формирователь M=LLT. Задайте 'ict' опция, чтобы использовать неполную факторизацию Холесского с пороговым отбрасыванием и задать диагональное значение сдвига 1e-6 избегать неположительных центров. Решите предобусловленную систему путем определения L и L' как вводит к symmlq.

setup = struct('type','ict','diagcomp',1e-6,'droptol',1e-14);
L = ichol(A,setup);
[x1,fl1,rr1,it1,rv1,rvcg1] = symmlq(A,b,tol,maxit,L,L');
fl1
fl1 = 0
rr1
rr1 = 2.7782e-15
it1
it1 = 3

Использование ichol предварительный формирователь значительно улучшает числовые свойства проблемы и symmlq может сходиться быстро. Выход rv1(1) norm(b), и выход rv1(end) norm(b-A*x1).

Можно следовать за прогрессом symmlq путем графического вывода относительных остаточных значений в каждой итерации. Постройте историю невязки метода сопряженных градиентов каждого решения с линией для заданного допуска.

semilogy(0:length(rvcg0)-1,rvcg0/norm(b),'-o')
hold on
semilogy(0:length(rvcg1)-1,rvcg1/norm(b),'-o')
yline(tol,'r--');
legend('No preconditioner','ICHOL preconditioner','Tolerance','Location','SouthEast')
xlabel('Iteration number')
ylabel('Relative residual')

Figure contains an axes. The axes contains 3 objects of type line, constantline. These objects represent No preconditioner, ICHOL preconditioner, Tolerance.

Исследуйте эффект предоставления symmlq с исходным предположением решения.

Создайте трехдиагональную разреженную матрицу. Используйте сумму каждой строки как вектор для правой стороны Ax=b так, чтобы ожидаемое решение для x вектор из единиц.

n = 900;
e = ones(n,1);
A = spdiags([e 2*e e],-1:1,n,n);
b = sum(A,2);

Используйте symmlq решить Ax=b дважды: одно время с исходным предположением по умолчанию, и одно время с хорошим исходным предположением решения. Используйте 200 итераций и допуск по умолчанию к обоим решениям. Задайте исходное предположение во втором решении как вектор со всеми элементами, равными 0.99.

maxit = 200;
x1 = symmlq(A,b,[],maxit);
symmlq converged at iteration 34 to a solution with relative residual 9.5e-07.
x0 = 0.99*e;
x2 = symmlq(A,b,[],maxit,[],[],x0);
symmlq converged at iteration 6 to a solution with relative residual 8.7e-07.

В этом случае, предоставляющем исходное предположение, включает symmlq сходиться более быстро.

Возвращение промежуточных результатов

Также можно использовать исходное предположение, чтобы получить промежуточные результаты путем вызова symmlq в цикле for. Каждый вызов решателя выполняет несколько итераций и хранит расчетное решение. Затем вы используете то решение в качестве начального вектора для следующего пакета итераций.

Например, этот код выполняет 100 итераций четыре раза и хранит вектор решения после каждой передачи в цикле for:

x0 = zeros(size(A,2),1);
tol = 1e-8;
maxit = 100;
for k = 1:4
    [x,flag,relres] = symmlq(A,b,tol,maxit,[],[],x0);
    X(:,k) = x;
    R(k) = relres;
    x0 = x;
end

X(:,k) вектор решения, вычисленный в итерации k из цикла for и R(k) относительная невязка того решения.

Решите линейную систему путем обеспечения symmlq с указателем на функцию, который вычисляет A*x вместо матрицы коэффициентов A.

Один из Уилкинсона тестирует матрицы, сгенерированные gallery 21 21 трехдиагональная матрица. Предварительно просмотрите матрицу.

A = gallery('wilk',21)
A = 21×21

    10     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     9     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     8     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1     7     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     1     6     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     1     5     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     1     4     1     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1     3     1     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     1     2     1     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     1     1     1     0     0     0     0     0     0     0     0     0     0
      ⋮

Матрица Уилкинсона имеет специальную структуру, таким образом, можно представлять операцию A*x с указателем на функцию. Когда A умножает вектор, большинством элементов в итоговом векторе являются нули. Ненулевые элементы в результате соответствуют ненулевым трехдиагональным элементам A. Кроме того, только основная диагональ имеет ненули, которые не равны 1.

Выражение Ax становится:

Ax=[1010001910001810017100161001510014100130001000110][x1x2x3x4x5x21]=[10x1+x2x1+9x2+x3x2+8x3+x4x19+9x20+x21x20+10x21].

Итоговый вектор может быть записан как сумма трех векторов:

Ax=[0+10x1+x2x1+9x2+x3x2+8x3+x4x19+9x20+x21x20+10x21+0]=[0x1x20]+[10x19x210x21]+[x2x210].

В MATLAB® запишите функцию, которая создает эти векторы и добавляет их вместе, таким образом давая значение A*x:

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end

(Эта функция сохранена как локальная функция в конце примера.)

Теперь решите линейную систему Ax=b путем обеспечения symmlq с указателем на функцию, который вычисляет A*x. Используйте допуск 1e-12 и 50 итераций.

b = ones(21,1);
tol = 1e-12;  
maxit = 50;
x1 = symmlq(@afun,b,tol,maxit)
symmlq converged at iteration 10 to a solution with relative residual 4.5e-16.
x1 = 21×1

    0.0910
    0.0899
    0.0999
    0.1109
    0.1241
    0.1443
    0.1544
    0.2383
    0.1309
    0.5000
      ⋮

Проверяйте тот afun(x1) дает вектор из единиц.

afun(x1)
ans = 21×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
      ⋮

Локальные функции

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end

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

свернуть все

Матрица коэффициентов в виде симметрической матрицы или указателя на функцию. Эта матрица является матрицей коэффициентов в линейной системе A*x = b. Обычно A большая разреженная матрица или указатель на функцию, который возвращает продукт большой разреженной матрицы и вектор-столбца. Можно использовать issymmetric подтвердить тот A issymmetric.

Определение A как указатель на функцию

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

Чтобы использовать указатель на функцию, используйте функциональную подпись function y = afun(x). Параметризация Функций объясняет, как предоставить дополнительные параметры функции afun, при необходимости. Вызов функции afun(x) должен возвратить значение A*x.

Типы данных: double | function_handle
Поддержка комплексного числа: Да

Правая сторона линейного уравнения в виде вектор-столбца. Длина b должно быть равно size(A,1).

Типы данных: double
Поддержка комплексного числа: Да

Допуск метода в виде положительной скалярной величины. Используйте этот вход для точности компромисса и времени выполнения в вычислении. symmlq должен соответствовать допуску в количестве позволенных итераций, чтобы быть успешным. Меньшее значение tol означает, что ответ должен быть более точным для вычисления, чтобы быть успешным.

Типы данных: double

Максимальное количество итераций в виде положительного скалярного целого числа. Увеличьте значение maxit позволить больше итераций для symmlq соответствовать допуску tol. Обычно меньшее значение tol средние значения больше итераций требуются, чтобы успешно завершать вычисление.

Матрицы перед формирователем в виде отдельных аргументов матриц или указателей на функцию. Можно задать матрицу перед формирователем M или его матричные факторы M = M1*M2 улучшить числовые аспекты линейной системы и облегчить для symmlq сходиться быстро. Можно использовать неполные матричные функции факторизации ilu и ichol сгенерировать матрицы перед формирователем. Также можно использовать equilibrate до факторизации, чтобы улучшить число обусловленности матрицы коэффициентов. Для получения дополнительной информации о предварительных формирователях смотрите Итерационные методы для Линейных систем.

symmlq обрабатывает незаданные предварительные формирователи как единичные матрицы.

Определение M как указатель на функцию

Можно опционально задать любой M, M1, или M2 как указатели на функцию вместо матриц. Указатель на функцию выполняет матрично-векторные операции вместо того, чтобы формировать целую матрицу перед формирователем, делая вычисление более эффективным.

Чтобы использовать указатель на функцию, используйте функциональную подпись function y = mfun(x). Параметризация Функций объясняет, как предоставить дополнительные параметры функции mfun, при необходимости. Вызов функции mfun(x) должен возвратить значение M\x или M2\(M1\x).

Типы данных: double | function_handle
Поддержка комплексного числа: Да

Исходное предположение в виде вектор-столбца с длиной равняется size(A,2). Если можно обеспечить symmlq с более разумным исходным предположением x0 чем нулевой вектор по умолчанию затем это может сохранить время вычисления и помочь алгоритму сходиться быстрее.

Типы данных: double
Поддержка комплексного числа: Да

Выходные аргументы

свернуть все

Решение для линейной системы, возвращенное как вектор-столбец. Этот выход дает приближенное решение линейной системы A*x = b. Если вычисление успешно (flag = 0), затем relres меньше чем или равно tol.

Каждый раз, когда вычисление не успешно (flag ~= 0), решение x возвращенный symmlq тот с минимальной нормой невязки, вычисленной по всем итерациям.

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

Флаговое значение

Сходимость

0

Успех — symmlq сходившийся к желаемому допуску tol в maxit итерации.

1

Отказ — symmlq выполненный с помощью итераций maxit итерации, но не сходились.

2

Отказ — матрица перед формирователем M или M = M1*M2 isIllConditioned.

3

Отказ — symmlq застоявшийся после того, как две последовательных итерации были тем же самым.

4

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

5

Отказ — матрица Перед формирователем M не симметричен положительный определенный.

Относительная остаточная ошибка, возвращенная как скаляр. Относительная остаточная ошибка является индикацией относительно как точный данный ответ x . symmlq отслеживает относительную невязку и невязку методов сопряженных градиентов в каждой итерации в процессе решения, и алгоритм сходится, когда любая невязка соответствует заданному допуску tol. relres выведите содержит значение невязки, которая сходилась, или относительная невязка или невязка методов сопряженных градиентов:

  • Относительная остаточная ошибка равна norm(b-A*x)/norm(b) и обычно невязка, которая соответствует допуску tol когда symmlq сходится. resvec выведите отслеживает историю этой невязки по всем итерациям.

  • Остаточная ошибка методов сопряженных градиентов равна norm(A'*A*x - A'*b). Эта невязка причины symmlq сходиться менее часто, чем относительная невязка. resveccg выведите отслеживает историю этой невязки по всем итерациям.

Типы данных: double

Номер итерации, возвращенный как скаляр. Этот выход указывает на номер итерации в который вычисленный ответ для x был вычислен.

Типы данных: double

Остаточная ошибка, возвращенная как вектор. Остаточная ошибка norm(b-A*x) показывает, как близко алгоритм к схождению для данного значения x. Число элементов в resvec равно количеству итераций. Можно исследовать содержимое resvec помочь решить, изменить ли значения tol или maxit.

Типы данных: double

Нормы невязки методов сопряженных градиентов, возвращенные как вектор. Число элементов в resveccg равно количеству итераций.

Типы данных: double

Больше о

свернуть все

Симметричный метод LQ

MINRES и методы SYMMLQ являются вариантами метода Lanczos, который подкрепляет метод сопряженных градиентов PCG. Как PCG, матрица коэффициентов все еще должна быть симметричной, но MINRES и SYMMLQ позволяют ему быть неопределенным (не, все собственные значения должны быть положительными). Это достигается путем предотвращения неявной LU-факторизации, обычно существующей в методе Lanczos, который подвержен отказам, когда с нулевыми центрами сталкиваются.

MINRES минимизирует невязку в 2-норме, в то время как SYMMLQ решает спроектированную систему с помощью факторизации LQ и сохраняет невязку ортогональной ко всем предыдущим единицам [1].

Советы

  • Сходимость большинства итерационных методов зависит от числа обусловленности матрицы коэффициентов, cond(A). Можно использовать equilibrate улучшить число обусловленности A, и самостоятельно это облегчает для большинства итеративных решателей сходиться. Однако использование equilibrate также приводит к лучшим качественным матрицам перед формирователем, когда вы впоследствии учитываете уравновешенный матричный B = R*P*A*C.

  • Можно использовать матричные функции переупорядочения такой как dissect и symrcm переставить строки и столбцы матрицы коэффициентов и минимизировать количество ненулей, когда матрица коэффициентов учтена, чтобы сгенерировать предварительный формирователь. Это может уменьшать память и время, требуемое впоследствии решить предобусловленную линейную систему.

Ссылки

[1] Барретт, R., М. Берри, Т. Ф. Чан, и др., Шаблоны для Решения Линейных систем: Базовые блоки для Итерационных методов, SIAM, Филадельфия, 1994.

[2] Пэйдж, C. C. и М. А. Сондерс, "Решение Разреженных Неопределенных Систем Линейных уравнений". SIAM J. Numer. Анальный., Vol.12, 1975, стр 617-629.

Расширенные возможности

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