Журнал вычислительной математики и математической физики, 2021, T. 61, № 5, стр. 787-799

1. ВВЕДЕНИЕ

Пусть $A\,:\lambda \mapsto {{\mathbb{C}}^{{n \times n}}}$ – аналитическое отображение на открытой области ${\text{\{ }}\lambda {\text{\} }} \subset \mathbb{C}$. Наша задача: найти $\lambda $ такое, что $f(\lambda ) = detA(\lambda ) = 0$.

Для вычисления нуля функции $f$ можно рассмотреть метод Ньютона:

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

Альтернативным подходом к вычислению производных, а значит, и ньютоновской коррекции, является алгоритмическое дифференцирование (см. [1]–[3]).

2. ВЫЧИСЛЕНИЕ ОПРЕДЕЛИТЕЛЕЙ И НЬЮТОНОВСКОЙ КОРРЕКЦИИ

При получении определителя обычно строится $LU$-разложение по методу Гаусса. Пусть $PA = LU$, где $P$ включает перестановки, связанные с частичным выбором ведущего элемента, $L$ – нижняя унитреугольная матрица и $U$ – верхняя треугольная матрица. Тогда

Когда $LU$-разложение записывается на месте матрицы $A$, текущее значение определителя на $k$-м шаге исключения умножается на ведущий элемент: $f: = f \times {{a}_{{kk}}}$. Это довольно быстро ведет к появлению очень больших или очень малых чисел. Поэтому лучше перейти к логарифмам: $logf: = logf + log({{a}_{{kk}}})$.

Заметим, что производная логарифма

является обратной величиной к ньютоновской коррекции. Таким образом, если производная логарифма вычисляется алгоритмическим дифференцированием то обратная величина ${\text{ffp}} = 1{\text{/logfp}} = f(\lambda ){\text{/}}f{\kern 1pt} {\text{'}}(\lambda )$ – в точности нужная нам ньютоновская коррекция. Это наблюдение используется в следующей программе (см. [2]):

function ffp=deta(A,Ap)

% DETA compute determinant of A and derivative

% Given A=A(lambda) and Ap=A'(lambda), DETA(A,Ap)

% computes Newton correction ffp=f/f' where f=det(A).

n=length(A); logfp=0;

for j=1:n

  [amax,kmax]= max(abs(A(j:n,j)));     % partial pivoting

  if amax == 0,ffp=0; return, end

  kmax=kmax+j-1;

  if kmax ~= j                 % interchange rows

    h=Ap(kmax,:); Ap(kmax,:)=Ap(j,:); Ap(j,:)=h;

    h=A(j,:); A(j,:)=A(kmax,:); A(kmax,:)=h;

  end

  logfp=logfp + Ap(j,j)/A(j,j);

  Ap(j+1:n,j)=(Ap(j+1:n,j)*A(j,j)-A(j+1:n,j)*Ap(j,j))/A(j,j)^2;

  A(j+1:n,j)=A(j+1:n,j)/A(j,j);

  Ap(j+1:n,j+1:n)=Ap(j+1:n,j+1:n) - Ap(j+1:n,j)*A(j,j+1:n)- …

             A(j+1:n,j)*Ap(j,j+1:n);

  A(j+1:n,j+1:n)=A(j+1:n,j+1:n) - A(j+1:n,j)*A(j,j+1:n);

end

ffp=1/logfp;

3. ПОДАВЛЕНИЕ ВМЕСТО ДЕФЛЯЦИИ

С помощью функции $deta$ мы можем вычислить решение уравнения $detA(\lambda ) = 0$ по методу Ньютона. Чтобы избежать перевычисления уже найденных нулей ${{\lambda }_{1}},\; \ldots ,\;{{\lambda }_{k}}$, мы подавим их, работая с функцией

где $p(\lambda ) = (\lambda - {{\lambda }_{1}})\; \cdots \;(\lambda - {{\lambda }_{k}})$. Тогда Производные для ${{f}_{k}}$ имеют вид (мы опускаем аргумент $\lambda $) Ньютонова коррекция ${{f}_{k}}{\text{/}}f_{k}^{'}$, выраженная через $f$ и $f{\kern 1pt} {\text{'}}$, приобретает вид В итоге получается итерация называемая итерацией Ньютона–Мехли (Newton–Maehly) (см. [4]).

В [2] мы привели расчеты для двух примеров масс с пружинами из [5] и для кубического примера из [1]. Вот $MATLAB$-программа для первого примера масс с пружинами:

n=50, tau=3, kappa=5,          % nonoverdamped

e=-ones(n-1,1);

C=(diag(e,-1)+diag(e,1)+3*eye(n)); K=kappa*C; C=tau*C;

lam=-0.5+0.1*i; lamb=[];        % start

for k=1:2*n

  ffp=1;

  while abs(ffp)>1e-14

    Qp=2*lam*eye(n)+C; Q=lam*(lam*eye(n)+ C)+K;

    ffp=deta(Q,Qp);

    s=sum(1./(lam-lamb(1:k-1)));

    lam=lam-ffp/(1-ffp*s);      % Newton step

  end

  lamb(k)=lam;

  lam=lam*(1+0.01*i);         % start for next eigenvalue

end

plot(lamb,'o')

Итерация для первого собственного значения начинается с выбора случайного комплексного числа, здесь ${{\lambda }_{0}} = - 0.5 + 0.1i$. В качестве начального значения для следующего собственного значения мы выбираем последнее из найденных и подавляем значение ${{\lambda }_{k}}$ с помощью малого возмущения: ${{\lambda }_{0}} = {{\lambda }_{k}}(1 + i{\text{/}}100)$.

Аналогично проводятся вычисления для перегруженной (overdamped) пружины и кубической проблемы собственных значений для $n = 50$. Найденные собственные значения изображены на фиг. 1.

Интересно посмотреть, сколько итераций нужно для каждого собственного значения. Столбцовые графики и средние числа итераций показаны на фиг. 2. Начальное значение для первого собственного значения, очевидно, выбрано не очень удачно, так как для сходимости требуется большое число итераций. Для кубической задачи большое число итераций возникает также для некоторых промежуточных собственных значений.

4. УЛУЧШЕНИЕ СХОДИМОСТИ

Причиной большого числа итерационных шагов при вычислении собственных значений в последних трех примерах является довольно плохая глобальная сходимость метода Ньютона. Локально метод Ньютона сходится к простому корню квадратично, обычно результат с машинной точностью получается за 3–4 итерации.

Пусть $f(z) = 0$ и ${{\lambda }_{k}}$ – приближение к $z$. Ньютонова итерация заменяет $f$ в окрестности ${{\lambda }_{k}}$ на линейную функцию $g$ такую, что $f({{\lambda }_{k}}) = g({{\lambda }_{k}})$, $f{\kern 1pt} {\text{'}}({{\lambda }_{k}}) = g{\text{'}}({{\lambda }_{k}})$. Таким образом, $g$ совпадает с отрезком ряда Тейлора $g(\lambda ) = f({{\lambda }_{k}}) + f{\kern 1pt} {\text{'}}({{\lambda }_{k}})(\lambda - {{\lambda }_{k}})$, а следующее приближение ${{\lambda }_{{k + 1}}}$ является нулем функции $g$.

Итерация Хэлли (Halley)

заменяет $f$ локально на гиперболическую функциютакую, что $f({{\lambda }_{k}}) = g({{\lambda }_{k}})$, $f{\kern 1pt} {\text{'}}({{\lambda }_{k}}) = g{\text{'}}({{\lambda }_{k}})$ и $f{\kern 1pt} {\text{''}}({{\lambda }_{k}}) = g{\text{''}}({{\lambda }_{k}})$, следующее приближение ${{\lambda }_{{k + 1}}}$ является нулем функции $g$. Итерация Хэлли – метод третьего порядка, и значит, он сходится к простому корню кубически (см. [6]). Можно ожидать, что гиперболическая аппроксимация для $f$ лучше с точки зрения глобальной сходимости.

5. РЕАЛИЗАЦИЯ ИТЕРАЦИИ ХЭЛЛИ

Нам нужна вторая производная определителя, более точно, нам нужно вычислять функцию

Заметим, что производная ньютоновской коррекции имеет вид Отсюда и чтобы получить $t(\lambda )$, нужно вычислять производную ньютоновской коррекции для нашей функции $deta$. Это делается алгоритмическим дифференцированием функции $deta$. Далее, функция $det2p$ получает на вход матрицы $A$, $A{\text{'}}$ и $A{\text{''}}$ и вычисляет ньютоновскую коррекцию $f{\text{/}}f{\kern 1pt} {\text{'}}$ и ее производную:

function [ffp,dffp] = det2p(A,Ap,App)

% DET2P computes Newton correction ffp = f/f'

%     and its derivative dffp = (f/f')'

n=length(A);

logfpp=0;                         % logfpp = log(f)''

logfp=0;                            % log(f)'

for k=1:n

  [amax,kmax]=max(abs(A(k:n,k)));    % partial pivoting

  if amax==0                       % matrix singular

   ffp=0; dffp=0;return

  end

  kmax=kmax+k-1;

  if kmax~=k                      % interchange rows

   h=App(k,:); App(k,:)=App(kmax,:); App(kmax,:)=h;

   h=Ap(k,:); Ap(k,:)=Ap(kmax,:); Ap(kmax,:)=h;

   h=A(k,:); A(k,:)=A(kmax,:); A(kmax,:)=h;

  end

  logfpp=logfpp+(A(k,k)*App(k,k)-Ap(k,k)^2)/A(k,k)^2;

  logfp=logfp+Ap(k,k)/A(k,k);

  App(k+1:n,k)=(A(k,k)*App(k+1:n,k)-Ap(k+1:n,k)*Ap(k,k))/A(k,k)^2-…

      (Ap(k+1:n,k)*Ap(k,k)/A(k,k)^2+ …

      A(k+1:n,k)*App(k,k)/A(k,k)^2-…

      2* A(k+1:n,k)*Ap(k,k)^2/A(k,k)^3);

  Ap(k+1:n,k)=Ap(k+1:n,k)/A(k,k)-A(k+1:n,k)*Ap(k,k)/A(k,k)^2;

  A(k+1:n,k)=A(k+1:n,k)/A(k,k);      % elimination step

  App(k+1:n,k+1:n)=App(k+1:n,k+1:n) -…

      (App(k+1:n,k)*A(k,k+1:n) + Ap(k+1:n,k)*Ap(k,k+1:n) ) - …

      (Ap(k+1:n,k)*Ap(k,k+1:n) + A(k+1:n,k)*App(k,k+1:n));

  Ap(k+1:n,k+1:n)=Ap(k+1:n,k+1:n) - Ap(k+1:n,k)*A(k,k+1:n)-…

      A(k+1:n,k)*Ap(k,k+1:n);

  A(k+1:n,k+1:n)=A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);

end

dffp=-logfpp/logfp^2; ffp=1/logfp;

6. ХЭЛЛИ–МЕХЛИ (HALLEY–MAEHLY)

Как и раньше, мы хотим подавить уже вычисленные собственные значения и снова рассматриваем

Теперь мы применяем итерации Хэлли к ${{f}_{k}}$: и записываем итерацию в терминах $f$, $f{\kern 1pt} {\text{'}}$ и $f{\kern 1pt} {\text{''}}$. Для ньютоновской коррекции ${{f}_{k}}{\text{/}}f_{k}^{'}$ используем уравнение (1). Для $f_{k}^{{''}}{\text{/}}f_{k}^{'}$ выражаем сначала Затем, деля на $f_{k}^{'}$, получаем и после умножения на ${{f}_{k}}{\text{/}}f_{k}^{'}$ находим Суммируя, реализуем итерацию Хэлли–Мехли, для этого необходимо следующее.

1. Вычислить ньютоновскую коррекцию для ${{f}_{k}}$:

2. Вычислить $t(\lambda )$ для ${{f}_{k}}$ в соответствии с уравнением (3).

3. Итерировать

Мы решаем три нелинейные проблемы собственных значений с помощью метода Хэлли и сравниваем результаты с итерациями Ньютона (см. табл. 1). Глобальная сходимость и в самом деле улучшилась, число итераций уменьшилось.

7. ЛАГЕР И ОСТРОВСКИЙ

Еще один метод третьего порядка для вычисления нулей многочлена – это метод Лагера. В качестве аппроксимации многочлена $f$ степени $n$ он использует многочлен $g(\lambda ) = a(\lambda - {{\lambda }_{1}}){{(\lambda - {{\lambda }_{2}})}^{{n - 1}}}$. Параметры $a$, ${{\lambda }_{1}}$ и ${{\lambda }_{2}}$ определяются таким образом, что $g$ интерполирует $f$ вместе с производными $f({{\lambda }_{k}}) = g({{\lambda }_{k}})$, $f{\kern 1pt} {\text{'}}({{\lambda }_{k}}) = g{\text{'}}({{\lambda }_{k}})$, $f{\kern 1pt} {\text{''}}({{\lambda }_{k}}) = g{\text{''}}({{\lambda }_{k}})$. Следующее приближение – это нуль для $g$, ближайший к ${{\lambda }_{k}}$:

Степень $n$ является параметром метода Лагера. Пусть в (4) $n \to \infty $. Тогда мы получаем итерацию которая называется итерацией Островского с квадратным корнем.

Заметим, что в итерационных методах Лагера и Островского так же, как и в методе Хэлли, нам нужно всего лишь два выражения:

Поскольку метод Лагера задумывался как метод вычисления нулей многочлена, можно ожидать, что он будет хорошо работать для трех наших примеров. Сравнение трех методов в табл. 2 показывает, что это действительно так.

8. МЕТОДЫ ТРЕТЬЕГО ПОРЯДКА

Методы Хэлли, Лагера и Островского реализуют специальные случаи следующей теоремы.

Теорема 1 (методы третьего порядка, см. [6]). Пусть $s$ – простой нуль функции $f$, а $G$ – любая функция такая, что

Тогда итерационный методсходится к $s$ по крайней мере кубически.

Примеры:

• формула Хэлли (Halley)

• формула Эйлера (Euler)

• квадратичная обратная интерполяция

• итерация Островского с квадратным корнем

• Лагер (Laguerre)

• семейство формул Хансен–Патрик (Hansen–Patrick) (см. [7])

Заметим, что для применения этих итераций нам нужно вычислять лишь ньютоновскую коррекцию и $t = ff{\kern 1pt} {\text{''/}}\mathop {f{\kern 1pt} '}\nolimits^2 $.

9. NLEVP – КОЛЛЕКЦИЯ НЕЛИНЕЙНЫХ ПРОБЛЕМ СОБСТВЕННЫХ ЗНАЧЕНИЙ

В замечательной коллекции NLEVP нелинейных проблем собственных значений есть примеры всех типов матриц (см. [8] и [9]). (http://www.maths.manchester.ac.uk/our-research/research-groups/numerical-analysis-and-scientific-computing/numerical-analysis/software/nlevp/)

Используя метод Лагера, мы решили две квадратичные задачи $sign1$ и $sign2$ (плотные матрицы, $n = 81$). Результаты приведены на фиг. 3. Задача $sign1$ имеет $2n = 162$ собственных значений на единичной окружности с двумя кластерами в точках $ \pm 1$. Из графиков следует, что сходимость к собственным значениям из этих двух кластеров медленная. Сходимость в задаче $sign2$ намного лучше, так как ее собственные значения лучше разделены. $MATLAB$-инструменты для замера времени $tic$, $toc$ показали, что на использованном автором ноутбуке задача $sign1$ решалась 63.92 с, а задача $sign2$ – 10.82 с.

10. НЕПОЛИНОМИАЛЬНАЯ ПРОБЛЕМА СОБСТВЕННЫХ ЗНАЧЕНИЙ

$TimeDelay$ представляет собой неполиномиальную нелинейную проблему собственных значений из NLEVP-коллекции с $3 \times 3$-матрицей $A(\lambda )$:

Об этой задаче пишут (см. [8]): “…характеристическое уравнение системы с запаздыванием по времени с единственной задержкой и постоянными коэффициентами. Задача имеет двойное непростое собственное значение $\lambda = 3\pi i$”.

Нелинейное  уравнение $detA(\lambda ) = 0$  имеет  бесконечно  много  решений. Используя итерации Островского,  мы можем,  например, вычислить первые 20 решений и получить график фиг.  4. На мнимой оси получаем упомянутые выше двойные собственные значения (${{\lambda }_{2}}$ и ${{\lambda }_{3}}$, фиг. 4) и также простое собственное значение ${{\lambda }_{4}} = 4.5\pi i$. Собственное значение ${{\lambda }_{1}} = 0.705244109106679 + 2.741466762205487i$ имеет положительную вещественную часть, остальные собственные значения имеют отрицательные вещественные части. Двойное собственное значение вычисляется с точностью, характерной для стандарта IEEE арифметики с плавающей точкой:

11. ИСКЛЮЧЕНИЕ ГАУССА ДЛЯ ЛЕНТОЧНЫХ МАТРИЦ

Многие задачи из NLEVP-коллекции имеют ленточные матрицы (например, $beamsensitivity$ с 7-ю или $pdde - stabiity$ с 32-мя диагоналями). Определители таких матриц разумно вычислять по специальному алгоритму. $MATLAB$-функция для исключения Гаусса для ленточных матриц с частичным выбором описана в [3]. Если $A$ имеет $q$ нижних и $p$ верхних диагоналей, то мы храним их как столбцы матрицы $B$. Для реализации частичного выбора добавляем $q$ нулевых столбцов к матрице $B$:

Адаптируя функцию $det2p$ к этой ленточной стуктуре, получаем функцию $det2pband$:

function [ffp,dffp]=det2pband(p,q,B,Bp,Bpp);

% DET2PBAND computes Newton-correction and derivative for a banded matrix

n=length(B); logfpp=0; logfp=0;

Bpp=[Bpp,zeros(n,q)];

Bp=[Bp,zeros(n,q)]; B=[B,zeros(n,q)];   % augment B with q columns

normb=norm(B,1);

for j=1:n

  maximum=0; kmax=j;                % search pivot

  for k=j:min(j+q,n)

    if abs(B(k,j-k+q+1))>maximum,

      kmax=k; maximum=abs(B(k,j-k+q+1));

    end

  end

  if maximum<1e-14*normb;             % only small pivots

    ffp=0; dffp=0; return               % consider det=0

  end

  if j~=kmax                       % interchange rows

    ind1=j-kmax+q+1:min(n,j+2*q+p-kmax+1);

    ind2=q+1:min(n,2*q+p+1);

    h=Bpp(kmax,ind1); Bpp(kmax,ind1)=Bpp(j,ind2); Bpp(j,ind2)=h;

    h=Bp(kmax,ind1); Bp(kmax,ind1)=Bp(j,ind2); Bp(j,ind2)=h;

    h=B(kmax,ind1); B(kmax,ind1)=B(j,ind2); B(j,ind2)=h;

  end

  logfpp=logfpp+(B(j,q+1)*Bpp(j,q+1)-Bp(j,q+1)^2)/B(j,q+1)^2;

  logfp=logfp+Bp(j,q+1)/B(j,q+1);

  for k=j+1:min(n,j+q)               % elimination step

    ind3=j-k+q+1;

    Bpp(k,ind3)=(Bpp(k,ind3)*B(j,q+1)-Bp(j,q+1)*Bp(k,ind3))/B(j,q+1)^2 …

            -(Bp(k,ind3)*Bp(j,q+1)+B(k,ind3)*Bpp(j,q+1))/B(j,q+1)^2 …

            +2*Bp(j,q+1)^2*B(k,ind3)/B(j,q+1)^3;

    Bp(k,ind3)=(B(j,q+1)*Bp(k,ind3)-B(k,ind3)*Bp(j,q+1))/B(j,q+1)^2;;

    B(k,ind3)=B(k,ind3)/B(j,q+1);

  end

  for k=j+1:min(n,j+q)

    for l=j+1:min(n,j+p+q)

      ind4=l-k+q+1; ind5=j-k+q+1; ind6=l-j+q+1;

      Bpp(k,ind4)=Bpp(k,ind4)-Bpp(k,ind5)*B(j,ind6)                -Bp(k,ind5)*Bp(j,ind6)…

              -Bp(k,ind5)*Bp(j,ind6)-B(k,ind5)*Bpp(j,ind6);

      Bp(k,ind4)=Bp(k,ind4)-Bp(k,ind5)*B(j,ind6)-B(k,ind5)*Bp(j,ind6);

      B(k,ind4)=B(k,ind4)-B(k,ind5)*B(j,ind6);

    end

  end

end

dffp=-logfpp/logfp^2; ffp=1/logfp;

Пример $beamsensitivity$ дает квадратичную проблему собственных значений с 7-диагональной матрицей. Для $n = 200$ вычисляются $400$ собственных значений по методу Лагера. Время измеряется $MATLAB$-функциями $tic$, $toc$. Если применить метод для плотной матрицы, то потребуется 83.85 с. Алгоритм с использованием ленточной структуры снижает время до 10.39 с.

12. ПРИМЕР СТРУНЫ В ВЯЗКОЙ СРЕДЕ

Хайям и др. пишут (см. [10]): “Стандартный подход к численному решению квадратичной задачи переводит $Q(\lambda ) = {{\lambda }^{2}}M + \lambda D + K$ в линейный полином $L(\lambda ) = \lambda X + Y$ с матрицей удвоенного размера с сохранением спектра. Уравнение $L(\lambda )z = 0$ обычно решается QZ-алгоритмом для задач умеренного размера и крыловскими методами для больших разреженных задач. Обычный выбор $L$ на практике – это первая сопровождающая форма вида

Когда $K$ и $M$ невырожденные, соответственно, два пучка: являются другими возможными линеаризациями”.

При использовании этих линеаризаций (см. [10]) результаты решения обобщенной проблемы собственных значений с помощью $MATLAB$-функции $eig$ довольно удручающие (фиг. 5).

Другие авторы (Fan, Lin, Van Dooren) показали в [11], как плохая обусловленность линеаризованных задач может быть исправлена масштабированием. В [10] объясняется, почему преобразованные системы без масштабирования являются настолько плохо обусловленными.

Ситуация напоминает мне старый пример Джима Уилкинсона, показывающий, что сведение проблемы собственных значений к задаче вычисления корней характеристического многочлена является не лучшим подходом к решению задачи, потому что изменяет ее обусловленность радикальным образом.

Однако при прямом решении уравнения $detA(\lambda ) = 0$ с помощью одного из наших методов, например, итераций Лагера, получаем корректные результаты без необходимости применять масштабирование (фиг. 6).

13. ВЫВОДЫ

Мы показали, как реализуются итерационные методы третьего порядка для решения $f(\lambda ) = detA(\lambda ) = 0$ с использованием автоматического (алгоритмического) дифференцирования. Эта техника дает точные производные, так как при вычислении определителей методом Гаусса используются лишь четыре арифметические операции.

Поскольку мы работаем непосредственно с исходной задачей, нет преобразований, которые могли бы изменить обусловленность задачи.

Кубическую сходимость можно получить с помощью многошаговых итераций, которые обходятся без вторых производных. Однако, используя только точечные итерации, мы получаем алгоритм $det2p$, в котором вычисляются ньютоновская коррекция и величина $t = ff{\kern 1pt} {\text{''/}}f{{{\kern 1pt} }^{{'2}}}$, содержащая нужные нам производные.

Вычисление вторых производных дорого. Для полной $n \times n$-матрицы нужно $ \sim {\kern 1pt} {{n}^{3}}$ операций на одну итерацию. Тем не менее, благодаря мощным процессорам наших компьютеров, мы можем решать нелинейные проблемы собственных значений умеренных размеров. Ситуация намного более благоприятна в случае ленточных матриц, для которых одна итерация требует лишь $ \sim {\kern 1pt} n$ операций.

Автор выражает благодарность Zhong-Zhi Bai и Yu-Mei Huang – организаторам Мемориального семинара в Ланжоу, посвященного Джину Голубу (Lanzhou, 2019). Приглашение участвовать в этой конференции дало особый стимул для разработки алгоритмов, обсуждаемых в этой статье.