Какие итеративные вспомогательные решатели, кроме бигга, можно использовать для решения уравнения Ньютона в MATLAB?

Question

Я пытаюсь проанализировать результаты для различных итерационных вспомогательных решателей системы уравнений Ньютона, используя переформулировку Ньютона-Фишера LCP (проблема линейной комплементарности). До сих пор я реализовал точный решатель - Gauss-Siedel и модифицированный Ньютоном метод, который использует bicg matlab как вспомогательный решатель уравнения J * h = -p (где J является jacobian, p - значение функции Фишера и h - это мой этап активации).

Часть коды, где я реализовал в BiCG subsolver и точный решатель:

if itt 
    % use iterative solver for newton eq 
    while ~all(fischer(x, A*x+b) == 0) & its < max_it 
     % compute the Jacobian 
     J = eval_jacobian(A, b, x); 
     % compute the value of the Fischer function 
     p = fischer(x, A*x + b); 
     % the natural merit function for convergence measure 
     residual(its) = .5*(p'*p); 
     % the newton eq, solve J*h = -p 
     h = bicg(J, -p, eps, s_its); 
     % update the solution vector 
     x = x + h; 
     % increment the iteration counter 
     its = its + 1; 
    end 
else 
    % the exact solver for newton equation 
    while ~all(fischer(x, A*x+b) == 0) & its < max_it 
     % compute the Jacobian 
     J = eval_jacobian(A, b, x); 
     % compute the value of the Fischer function 
     p = fischer(x, A*x + b); 
     % the natural merit function for convergence measure 
     residual(its) = .5*(p'*p); 
     % the newton eq, solve J*h = -p 
     h = - J/p; 
     % update the solution vector 
     x = x + h; 
     % increment the iteration counter 
     its = its + 1; 
    end 

Итак, мой вопрос будет то, что другой итерационный суб-решателей вы бы использовали? Я не возражаю против их реализации, если для них нет функции в matlab. Я понимаю, что это скорее теоретический вопрос. Спасибо.

Answer 1

В дополнение к bicg существует целый ряд других итерационных методов для общих несимметричных систем, может иногда предложение лучше конвергенции, такие как bicgstab и gmres. Оба эти метода аналогичны по духу bicg, но связаны с различными стратегиями обновления в подпространстве крылова. MATLAB имеет реализации для обоих этих методов, поэтому вы можете проверить их на свою проблему.

В общем, при использовании итерационных решателей крылова вы часто можете добиться лучшей производительности, включив предварительное кондиционирование. Это довольно привлекательный предмет - если вы еще не знакомы, here - это хорошее место для начала. Я бы посоветовал начать с простого предварительного кондиционирования jacobi diag(J) и перейти оттуда.

Итерационные решатели также могут получить дополнительную эффективность на основе так называемого «без матричного» подхода - решателю на самом деле нужно только вычислить вектор-векторный продукт J*p, ему не нужна полная матрица J. Может быть возможно написать эффективную подпрограмму, которая вычисляет J*p без форматирования J, но это будет зависеть от деталей конкретной проблемы.

Последнее соображение - сам программирующий язык - MATLAB использует высоко оптимизированный библиотечный код для его прямых решателей J\p, тогда как методы крылова реализованы как m-код, что обычно приводит к штрафу за производительность.

Надеюсь, это поможет.