Алгоритм решения задачи на собственные значения

[Аурелиано Буэндиа]

Нужен эффективный алгоритм решения обобщенной задачи на собственные значения
${\bf A}x=\lambda{\bf B}x$,

где ${\bf A},{\bf B}$ - эрмитовые матрицы $n\times n$, а $x$ - вектор-столбец, ${\lambda}$ - собственное значение задачи
Никто не сталкивался с такими задачами?

[maxal]

Если матрица B обратима, то задача сводится к нахождению собственных чисел матрицы B^{-1} A:
(B^{-1} A) x = \lambda x

Если матрица A обратима, то задача сводится к нахождению собственных чисел матрицы A^{-1} B:
(A^{-1} B) x = \lambda^{-1} x

Единственный нетривиальный случай, когда обе матрицы A и B необратимы.

[Yuri Gendelman]

В этой книге: Уилкинсон, Райнш. Справочник алгоритмов на языке АЛГОЛ: Линейная алгебра. М., "Машиностроение", 1976, приведен алгоритм для симметричных вещественных матриц (II.10. стр.267-276):

1. Треугольное разложение по Холецкому: B = L L^{T}
2. Сведение к стандартной симмертичной задаче на С.З.
(L^{-1} A L^{-T}) (L^{T} x) = \lambda (L^{T} x)
Получаем С.З. и "повернутые С.В." (y)
3. Решение системы L^{T} x = y

Можно попробовать заменить п.1 на что-нибудь типа B = U^{-1} Q U, где U - унитарная, а Q - диагональная (и вещественная), и получить модификацию алгоритма для эрмитовых матриц:
2. (U A U^{-1}) (U x) = (\lambda Q) (U x)
3. Решить U x = y и \lambda_{i} = \lambda^{mod}_{i} / Q_{i}_{i}

[Аурелиано Буэндиа]

maxal писал(а):

Если матрица B обратима, то задача сводится к нахождению собственных чисел матрицы B^{-1} A:
(B^{-1} A) x = \lambda x

Произведение эрмитовых матриц может быть неэрмитовым. Так и получается. В этом случае нужен уже алгоритм для неэрмитовых матриц.

У меня эрмитовые комплексные матрицы

[maxal]

Аурелиано Буэндиа писал(а):

В этом случае нужен уже алгоритм для неэрмитовых матриц.

А что в этом такого плохого?

[Аурелиано Буэндиа]

maxal писал(а):

А что в этом такого плохого?

Ничего. Только нужен стабильно работающий алгоритм.

[Аурелиано Буэндиа]

Yuri Gendelman писал(а):

В этой книге: Уилкинсон, Райнш. Справочник алгоритмов на языке АЛГОЛ: Линейная алгебра. М., "Машиностроение", 1976, приведен алгоритм для симметричных вещественных матриц (II.10. стр.267-276):

1. Треугольное разложение по Холецкому: B = L L^{T}
2. Сведение к стандартной симмертичной задаче на С.З.
(L^{-1} A L^{-T}) (L^{T} x) = \lambda (L^{T} x)
Получаем С.З. и "повернутые С.В." (y)
3. Решение системы L^{T} x = y

Можно попробовать заменить п.1 на что-нибудь типа B = U^{-1} Q U, где U - унитарная, а Q - диагональная (и вещественная), и получить модификацию алгоритма для эрмитовых матриц:
2. (U A U^{-1}) (U x) = (\lambda Q) (U x)
3. Решить U x = y и \lambda_{i} = \lambda^{mod}_{i} / Q_{i}_{i}

Yuri Gendelman, я забыл вас поблагодарить за попытку помочь. Вообще идея с разложением по Холецкому мне понравилась. Но в моем случае в задаче ${\bf A}x=\lambda {\bf B}x$ матрицы ${\bf A},{\bf B}$ не являются положительно определенными. и насколько я знаю разложение по Холецкому для этого случая не существует. А что касается 3-го пункта, то допустим я сделал разложение B = U^{-1} Q U, но откуда следует что \lambda_{i} = \lambda^{mod}_{i} / Q_{i}_{i}?
Прокомментируйте пожалуйста.

[Yuri Gendelman]

Всегда пожалуйста.

А насчет Холецкого - это недоразумение.
1. Я сначала описал опубликованный алгоритм для симметричных вещественных (еще уточню - положительно определенных) матриц, где используется разложение Холецкого.
2. Затем я предложил модификацию для эрмитовых матриц. При этом разложение Холецкого разумеется не подходит. Нужно разложение вида $B = U^{-1} Q U$, где U - унитарная, а Q - диагональная (и вещественная) матрица, то есть $Q_{i,j} = \delta_{i,j} Q_{i,i}$ . Такое разложение существует и используется в алгоритмах для стандартной задачи на с.з.с.в. для эрмитовых матриц.
3. При решении $(U A U^{-1}) (U x) = (\lambda Q) (U x)$ произведение $(\lambda Q) == \lambda^{mod}$ - это произведение двух диагональных матриц , где на диагонали - $\lambda_{i} Q_{i,i}$. Отсюда - интересующая Вас формула.

Если мои объяснения показались не вполне ясными - welcome. Я этим алгоритмом не пользовался и правильность не доказал.
Кстати, Вам наверное известна библиотека Библиотека Численного Анализа НИВЦ МГУ (БЧА НИВЦ МГУ) http://www.srcc.msu.su/num_anal/lib_na/cat/cat0.htm ?

[Аурелиано Буэндиа]

Yuri Gendelman писал(а):

При решении $(U A U^{-1}) (U x) = (\lambda Q) (U x)$ произведение $(\lambda Q) == \lambda^{mod}$ - это произведение двух диагональных матриц , где на диагонали - $\lambda_{i} Q_{i,i}$. Отсюда - интересующая Вас формула.

У меня большие сомнения по поводу правильности такого подхода. Есть контрпример с
матрицами
\left( \begin{array}{cc} 0 & 1 \\ 1 & 0 \\ \end{array} \right) x = \lambda \left( \begin{array}{cc} 1 & 0 \\ 0 & -4 \\ \end{array} \right) x

Yuri Gendelman писал(а):

Кстати, Вам наверное известна библиотека Библиотека Численного Анализа НИВЦ МГУ (БЧА НИВЦ МГУ) http://www.srcc.msu.su/num_anal/lib_na/cat/cat0.htm ?

Теперь известна. Спасибо.

[Yuri Gendelman]

Аурелиано Буэндиа писал(а):

У меня большие сомнения по поводу правильности такого подхода.

Да, Вы правы, я ошибся.

[Comga]

Аурелиано Буэндиа писал(а):

Нужен эффективный алгоритм решения обобщенной задачи на собственные значения
${\bf A}x=\lambda{\bf B}x$,

где ${\bf A},{\bf B}$ - эрмитовые матрицы $n\times n$, а $x$ - вектор-столбец, ${\lambda}$ - собственное значение задачи
Никто не сталкивался с такими задачами?

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

[Аурелиано Буэндиа]

Comga писал(а):

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

Меня пока интересует метод для $n = 100$

[Comga]

Аурелиано Буэндиа писал(а):

Меня пока интересует метод для $n = 100$

Для малых размерностей посоветую начать с Eispack, Lapack , для больших с Y.Saad. Numerical Methods for Large EigenValue Problems. Заодно и нам расскажите.

[photon]

Возможно, чем-то поможет это:
http://www.netlib.org/lapack/lug/lapack_lug.html

А еще я нашел, я не знаю, что это такое, но MatLAB именно для случая эрмитовых матриц по умолчанию использует "Cholesky factorization" - возможно, стоит поискать по этому ключу.

[Аурелиано Буэндиа]

Comga, photon , спасибо.