Собственные числа матрицы

[Joyce]

Здравствуйте! Столкнулся со следующим вопросом: есть матрица(ы) с рациональными числами, матрица(ы) самосопряжённая, т.е. собственные числа действительные. Хотелось бы их посдчитать приблежённо. Но, при этом понимать какая там погрешность, вот, сколько точных цифр после запятой в десятичном виде.

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

И было бы замечательно, если бы кто-нибудь указал на обзор о таких вещах, как люди это всё дело считают. Сам я с такой задачей, о приближённом вычислении, впервые столкнулся, так что сильно не пинайте.

[zykov]

Если матрица небольшая, то составьте характеристический многочлен (например в wxMaxima легко сделать). Его коэффициенты будут тоже рациональными. А если привести к общему знаменателю, то будет многочлен с целыми коэффициентами.
Дальше можете аппроксимировать его корни рациональными числами с любой точностью (например методом Ньютона).

[мат-ламер]

Joyce в сообщении #1537062 писал(а):

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

Я пользуюсь MatLab-ом и Maple. Если матрица сильно большого размера, то возможно MatLab предпочтительней. В документации они иногда приводят ссылку на статью с алгоритмом. Но в это лучше не лезть. Точность вычислений вообще сложно определить. Вообще есть понятие обусловленности специально для задачи нахождения собственных значений. Но эту обусловленность непонятно как считать. Можно попробовать несколько раз запустить счёт, меняя при этом исходные данные. И смотреть при этом, как меняется ответ.

[Joyce]

zykov в сообщении #1537064 писал(а):

Дальше можете аппроксимировать его корни рациональными числами с любой точностью (например методом Ньютона).

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

И не очень понятно для меня, матрица $51\times51$ -- большая или нет.

-- 30.10.2021, 18:45 --

Спасибо всем за ответы.

мат-ламер в сообщении #1537067 писал(а):

Точность вычислений вообще сложно определить. Вообще есть понятие обусловленности специально для задачи нахождения собственных значений.

Т.е. точность, как я понял, определить сложно. Жаль. Меня этот вопрос интересует ещё и в том разрезе, что видел когда-то статьи, где к аналитическим задачам применяют численные методы и что-то приближённо считают. Как авторы обосновывают свои приближённые ответы для меня интересно. Я так предполагаю, что это очень обширная область. Но лезть туда, желания нет.

[zykov]

Ну лучше 10x10. Но и 50x50 вполне в рамках.
Один корень можно и вручную там же в wxMaxima. Боюсь, 50 корней вручную - слишком.

Тут многое зависит от того, что надо.
Эти рациональные значения в матрице - это точные значения или они сами являются приближениями?
Вообще wxMaxima и собственные значения, и вектора сама считает. Может и точно посчитать, но это для маленьких матриц, где корни многочлена можно точно выразить.

Если значения были приближениями, то вопрос с какой точностью?
Если для их представления достаточно double float, то да, в Matlab/Octave проще просто численно посчитать.

[Joyce]

zykov спасибо, посмотрю эту wxMaxima!

zykov в сообщении #1537069 писал(а):

Эти рациональные значения в матрице - это точные значения или они сами являются приближениями?

Значения матрицы являются точными. У меня явный вид есть.

[zykov]

Joyce в сообщении #1537074 писал(а):

Значения матрицы являются точными.

Хорошо, тогда и многочлен будет точным. Тогда проще проанализировать точность приближений его корней.
Вот пример кода для Maxima:

Используется синтаксис Text

(%i19)  m: matrix(
         [1/2,6/7,0,3,17], 
         [6/7,1/3,7/8,0,3], 
         [0,7/8,1/4,8/9,0], 
         [3,0,8/9,1/5,9/10], 
         [17,3,0,9/10,1/6]
        );
(m)     matrix(
                [1/2,   6/7,    0,      3,      17],
                [6/7,   1/3,    7/8,    0,      3],
                [0,     7/8,    1/4,    8/9,    0],
                [3,     0,      8/9,    1/5,    9/10],
                [17,    3,      0,      9/10,   1/6]
        )
(%i20)  f:expand(charpoly(m, x));
(f)     -x^5+(29*x^4)/20+(1964146249*x^3)/6350400-(412926571*x^2)/6350400-(63388008323*x)/114307200+24274979531/91445760
(%i21)  bfallroots(f);
(%o21)  [x=5.283621589070233b-1,x=-1.451760794091658b0,x=1.124562731604294b0,x=-1.692922925739396b1,x=1.81780651609743b1]
 

Многочлен будет точный (с рациональными коэффициентами).
'bfallroots' ищет корни многочлена приближенно, но использует bigfloat, для которых можно задать произвольную точность (подробнее есть в документации Maxima, переменная 'fpprec' задаёт точность).

-- 30.10.2021, 20:55 --

Собственно, все цифры десятичной дроби должны быть точными.
Задать количество этих цифр (например 40 цифр) можно так:

Используется синтаксис Text

fpprec:40

[мат-ламер]

Joyce в сообщении #1537074 писал(а):

Значения матрицы являются точными. У меня явный вид есть.

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

[Евгений Машеров]

Если есть и собственные вектора, и собственные значения - восстановить исходную матрицу и сравнить с изначальной.

[zykov]

Попробовал этот подход в wxMaxima, к сожалению больше 10x10 не тянет.
Придется численно спектр искать. Например в Matlab или его бесплатном аналоге Octave.
Как тут уже писали, из алгоритма трудно точность понять. Проще результат проанализировать, насколько он точен.

[Евгений Машеров]

Вопрос о точности оценок собственных значений рассматривался давно, в классической книге Уилкинсона "Алгебраическая проблема собственных значений", а применительно к симметричным - Парлетт, "Симметричная проблема собственных значений" (обе книги переведены на русский, и достаточно давно, но я, ограничиваясь прикладным вопросом - насколько велика вычислительная погрешность сравнительно с точностью исходных данных, ограничился и тамошники результатами, а так можно поискать и более новые результаты).
Практический же мой вывод состоял в том, что существующие методы дают достаточно высокую точность, причём такой неприятности, как, скажем, при решении СЛАУ с матрицей, близкой к вырожденной, не бывает, не возникает "усиления ошибки".
Если Ваш метод даёт и собственные значения $\lambda$, и собственные вектора $x$, проще всего подставить $Ax=\lambda x$, или скорее $y=Ax$, а затем посчитать норму игрека и насколько она отличается от соответствующей лямбды.
Если с.в. недоступны, можно посчитать $A-\lambda_i I$, а затем ортогонализовать одну из строк к прочим. При этом в точном случае должен быть ноль, а если не ноль, то, значит, норма ортогонализированной строки и покажет ошибку.
Через характеристический многочлен лично я бы поостерёгся, расчёт "определителя по определению" это факториальная сложность, и где-то после $n=10$ уже бессмысленная трата машинного времени. В своё время, правда, были созданы многочисленные способы нахождения коэффициентов характеристического многочлена с полиномиальной сложностьб (Фадеев и Фадеева), но возникает вопрос о нахождении ошибки этих методов.

[svv]

Есть такая книга (можно скачать на LibGen):
Годунов, Антонов, Кирилюк, Костин. Гарантированная точность решения систем линейных уравнений в евклидовых пространствах. Новосибирск, 1988.

Что сделали авторы. Они предложили такой метод решения СЛАУ (метод ортогональных преобразований), для которого возможно строго оценить его точность. Ну, и оценили. Сделали это честно: одна из глав называется «Особенности машинных вычислений».

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

Книга вообще непростая, у меня бы ушёл год на её изучение от начала до конца. Но хороший алгоритм решения СЛАУ (так сказать, минимально портящий матрицу в процессе решения) я оттуда позаимствовал, реализовал на C++ и потом много где использовал.

[zykov]

zykov в сообщении #1537099 писал(а):

Попробовал этот подход в wxMaxima, к сожалению больше 10x10 не тянет.

Немного повозился и удалось заставить работать.
Функции 'charpoly' и 'determinant' плохо работают, если матрица не маленькая (больше 10x10). Видимо сделаны "в лоб" без оптимизации.
Зато вышло найти детерминант через функцию 'lu_factor' и построить точный характеристический многочлен.
Вот код для Maxima:

Используется синтаксис Text

n:51;
m: genmatrix(lambda([i,j], 51-abs(i-j)), n, n);
res:[];
for x:1 thru n+1 do
  (r:lu_factor(m-x*ident(n)),
  p:1,
  for i:1 thru n-1 do
    for i1:i+1 thru n do (p:p*(if r[2][i1]<r[2][i] then -1 else 1)),
  for i:1 thru n do (p:p*r[1][r[2][i]][i]),
  res:endcons(p, res));
res;
b: genmatrix(lambda([i,j], res[i]), n+1, 1);
am: genmatrix(lambda([i,j], i^(j-1)), n+1, n+1);
pol1: invert(am).b;
b1: genmatrix(lambda([i,j], x^(i-1)), n+1, 1);
f: transpose(pol1).b1;

rt:realroots(f, 1e-25);
 

Его надо в файл записать и вызывать как

Используется синтаксис Text

batch("полное_имя_файла");

Во второй строке генерируется некоторая матрица 51x51. Это надо заменить на нужную матрицу.
В 'f' будет сам характеристический многочлен.
В конце 'realroots' ищет приближения его корней в виде рациональных дробей с точностью указанной во втором параметре (можно хоть 1e-100 указать).
При желании потом эти рациональные дроби можно в big float преобразовать как 'bfloat(rt);', предварительно выставив 'fpprec:25'.

[Евгений Машеров]

Всё замечательно, но считать через характеристический многочлен - не самый точный способ...

[vpb]

zykov в сообщении #1537064 писал(а):

Если матрица небольшая, то составьте характеристический многочлен (например в wxMaxima легко сделать). Его коэффициенты будут тоже рациональными. А если привести к общему знаменателю, то будет многочлен с целыми коэффициентами.
Дальше можете аппроксимировать его корни рациональными числами с любой точностью (например методом Ньютона).

Не в обиду будь сказано, Ваш совет совершенно ошибочен.

Перейдем, однако, к делу.
Joyce
Уточним, чего вам надо. Верно ли, что (а) вы решаете данный вопрос для, фактически, однократного или близкого к тому использования ? Ваша программа не на продажу, и не для того, чтобы применять ее автоматически миллиарды раз ? (В последнем случае имела бы значение общая производительность программы; но для однократного применения в принципе всё равно, будут с.з. матрицы размера 50 на 50 считаться одну миллисекунду илм одну минуту). Просто есть один конкретный самосопряженный оператор, размерности около 50, и вам надо более-менее точно определить его спектр ?

Если так, задача решается довольно просто. Дело в том, что спектр симметрической матрицы очень устойчив к небольшим возмущениям элементов. Это можно уточнить по-разному, например в виде т.наз. теоремы Виландта-Гофмана. См., например, книгу Уилкинсона. И к возможным ошибкам округления при реализации какого-либо разумного алгоритма тоже устойчив. Короче, следствие такое: можете брать любой пакет, в котором зашито вычисление собственных значений симметрической матрицы, и результат будет очень точен. То есть, в вашем случае, думаю, если матрица размера 50, и все вычисления делаются с точностью double, т.е. с относительной точностью около $10^{-16}$, то ошибка в любом
собственном значении будет не более $10^{-13}M$, где $M$ --- наибольшее по модулю точное собственное значение. Это можно доказать (я в уме прикинул...), а фактически ошибка скорее всего еще меньше будет. А если хотите точность порядка последней разрядной единицы в наибольшем с.з., то можно делать вычисления, скажем, с точностью $10^{-22}$ или с большей, а потом найденные собственные значения округлить до double.

Для несимметрических же матриц ситуация совершенно другая, их с.з. считать гораздо труднее. Насчет них я не в курсе.

Впрочем, и искать какой-то пакет необязательно. Я, например, года три назад реализовал сам на Сях метод вращений Якоби, для своих собственных надобностей, и был вполне доволен результатом. И это немного времени заняло, пару дней (если не один вообще). Притом сам я отнюдь не программист.

Теперь насчет чего почитать по численной линейной алгебре. Гм. Курочка по зернышку клюет, поэтому я перечислю просто те книжки, из которых прочел более-менее много (из каких двадцать страниц, а из каких и двести). (Может тоже что для себя наклюете... ).

1) Бахвалов, Численные методы, глава про численные методы алгебры

2) Воеводин, Вычислительные основы линейной алгебры

3) Беклемишев, Дополнительные главы линейной алгебры (или как она там называется...)

4) Wilkinson, Rounding errors in algebraic processes

5) Форсайт и Молер, Численное решение систем линейных алгебраических уравнений

6) Overton, Numerical computing with IEEE floating point arithmetic

7) Higham, Accuracy and stability of numerical algorithms

8) Уилкинсон, Алгебраическая проблема собственных значений

9) Trefethen, Bau, Numerical linear algebra

10) Деммель, Вычислительная линейная алгебра

11) Годунов с соавторами, упомянута выше в теме

12) Парлетт, упомянута там же

Еще есть книги (не читал практически)

13) Голуб, ван Лоан, Матричные вычисления

14) Уилкинсон, Райнш, Справочник алгоритмов на языке Алгол

15) Фаддеев и Фаддева, Вычислительная линейная алгебра (в настоящее время представляет
лишь исторический интерес),

16) Малышев, Введение в вычислительную линейную алгебру.

Из них, наиболее важны были Уилкинсон (малый, номер 4), Форсайт и Молер, Овертон (это
просто мастхев, как нынче говорят; собственно, она не про линейную алгебру, а про компьютерную
арифметику), и Трефетен-Бау.

-- 01.11.2021, 05:38 --

zykov в сообщении #1537113 писал(а):

m: genmatrix(lambda([i,j], 51-abs(i-j)), n, n);

не понял. Я с синтаксисом Максимы не знаком, системы компьютерной алгебры не пользую. Вы взяли матрицу с элементами $\lambda_{ij}=51-|i-j|$, да ? Но, такая матрица, если я не путаю (давненько не брал я в руки шашек!), более-менее тривиально решается руками, и в качестве тестовой не подходит. Попробуйте что-нибудь другое ... ну , например, обратные величины чисел от $1$ до $51^2$, расположенные в естественном порядке (как в книжке). И посмотрите, не издохнет ли эта ваша максима.