Численные методы вычисления собственных значений матриц

Введение

Численные методы вычисления собственных значений матрицы
В численном анализе одной из важнейших задач является разработка эффективных и устойчивых алгоритмов для нахождения собственных значений матрицы. Эти алгоритмы для собственных значений могут также находить собственные векторы.

Нормальные, гермитовые и реально-симметричные матрицы

Присоединенная M^(*) комплексной матрицы M является транспонированием сопряжения M: квадратная матрица A называется нормальной, если она коммутирует со своей присоединенной: A^(*)A = AA^(*). Она называется эрмитовой, если она равна своей присоединенной: A^(*) = A. Все эрмитовы матрицы являются нормальными. Если A имеет только вещественные элементы, то присоединенная – это просто транспонирование, а A является эрмитовой тогда и только тогда, когда она симметрична. При применении к столбцовым векторам, присоединенная может быть использована для определения канонического внутреннего произведения на C^n: <w, v> = w^(*)v. Нормальные, эрмитовы и вещественно симметричные матрицы обладают несколькими полезными свойствами: каждый обобщенный собственный вектор нормальной матрицы является обычным собственным вектором. Любая нормальная матрица подобна диагональной матрице, поскольку ее жорданова нормальная форма диагональна. Собственные векторы различных собственных значений нормальной матрицы ортогональны. Нулевое пространство и образ (или пространство столбцов) нормальной матрицы ортогональны друг другу. Для любой нормальной матрицы A, C^n имеет ортонормированный базис, состоящий из собственных векторов A. Соответствующая матрица собственных векторов унитарна. Собственные значения эрмитовой матрицы являются вещественными, поскольку для ненулевого собственного вектора v это справедливо. Если A вещественная, то существует ортонормированный базис для R^n, состоящий из собственных векторов A, тогда и только тогда, когда A симметрична. Реальная или комплексная матрица может иметь все вещественные собственные значения, не являясь при этом эрмитовой. Например, вещественная треугольная матрица имеет свои собственные значения на главной диагонали, но в общем случае не является симметричной.

Номер условия

Любую задачу численного вычисления можно рассматривать как вычисление некоторой функции f для заданного x. Условное число κ(f, x) задачи – это отношение относительной погрешности в выходных данных функции к относительной погрешности во входных данных, и оно зависит как от функции, так и от входных данных. Условное число описывает, как погрешность нарастает в процессе вычисления. Его десятичный логарифм показывает, на сколько меньше значащих цифр в результате, чем во входных данных. Условное число представляет собой наилучший случай. Оно отражает присущую задаче неустойчивость, независимо от способа её решения. Ни один алгоритм не может дать более точных результатов, чем указано условным числом, за исключением случайности. Однако, плохо разработанный алгоритм может привести к значительно худшим результатам. Например, как упоминается ниже, задача нахождения собственных значений для нормальных матриц всегда хорошо обусловлена. Однако, задача поиска корней многочлена может быть очень плохо обусловлена. Следовательно, алгоритмы собственных значений, работающие путем нахождения корней характеристического многочлена, могут быть плохо обусловлены, даже если сама задача хорошо обусловлена. Для задачи решения линейного уравнения 1=A'v' = 'b', где A обратима, матричное условное число 1=κ(A^(−1), 'b') задается выражением , где – операторная норма, подчиненная обычной евклидовой норме на 'C'^(n). Поскольку это число не зависит от 'b' и одинаково для A и A^(−1), его обычно называют просто условным числом κ(A) матрицы A. Это значение κ(A) также является абсолютным значением отношения наибольшего собственного значения A к наименьшему. Если A унитарна, то , следовательно 1=κ(A) = 1. Для общих матриц операторную норму часто трудно вычислить. По этой причине для оценки условного числа обычно используются другие матричные нормы. Для задачи собственных значений Бауэр и Фик доказали, что если λ является собственным значением диагонализуемой n × n матрицы A с матрицей собственных векторов V, то абсолютная погрешность в вычислении λ ограничена произведением κ(V) и абсолютной погрешности в A. В результате, условное число для нахождения λ равно . Если A нормальна, то V унитарна, и 1=κ(λ, A) = 1. Таким образом, задача собственных значений для всех нормальных матриц хорошо обусловлена. Условное число для задачи нахождения собственного пространства нормальной матрицы A, соответствующего собственному значению λ, было показано обратно пропорциональным минимальному расстоянию между λ и другими различными собственными значениями A. В частности, задача собственного пространства для нормальных матриц хорошо обусловлена для изолированных собственных значений. Когда собственные значения не изолированы, лучшее, на что можно надеяться, – это определить линейную оболочку всех собственных векторов близлежащих собственных значений.

Any problem of numeric calculation can be viewed as the evaluation of some function f for some input x. The condition number κ(f, x) of the problem is the ratio of the relative error in the function's output to the relative error in the input, and varies with both the function and the input. The condition number describes how error grows during the calculation. Its base 10 logarithm tells how many fewer digits of accuracy exist in the result than existed in the input. The condition number is a best case scenario. It reflects the instability built into the problem, regardless of how it is solved. No algorithm can ever produce more accurate results than indicated by the condition number, except by chance. However, a poorly designed algorithm may produce significantly worse results. For example, as mentioned below, the problem of finding eigenvalues for normal matrices is always well conditioned. However, the problem of finding the roots of a polynomial can be very ill conditioned. Thus eigenvalue algorithms that work by finding the roots of the characteristic polynomial can be ill conditioned even when the problem is not. For the problem of solving the linear equation 1=A'v' = 'b' where A is invertible, the matrix condition number 1=κ(A^(−1), 'b') is given by , where is the operator norm subordinate to the normal Euclidean norm on 'C'^(n). Since this number is independent of 'b' and is the same for A and A^(−1), it is usually just called the condition number κ(A) of the matrix A. This value κ(A) is also the absolute value of the ratio of the largest eigenvalue of A to its smallest. If A is unitary, then , so 1=κ(A) = 1. For general matrices, the operator norm is often difficult to calculate. For this reason, other matrix norms are commonly used to estimate the condition number. For the eigenvalue problem, Bauer and Fike proved that if λ is an eigenvalue for a diagonalizable n × n matrix A with eigenvector matrix V, then the absolute error in calculating λ is bounded by the product of κ(V) and the absolute error in A. As a result, the condition number for finding λ is If A is normal, then V is unitary, and 1=κ(λ, A) = 1. Thus the eigenvalue problem for all normal matrices is well conditioned. The condition number for the problem of finding the eigenspace of a normal matrix A corresponding to an eigenvalue λ has been shown to be inversely proportional to the minimum distance between λ and the other distinct eigenvalues of A. In particular, the eigenspace problem for normal matrices is well conditioned for isolated eigenvalues. When eigenvalues are not isolated, the best that can be hoped for is to identify the span of all eigenvectors of nearby eigenvalues.

Итеративные алгоритмы

Итеративные алгоритмы решают задачу собственных значений, генерируя последовательности, сходящиеся к собственным значениям. Некоторые алгоритмы также генерируют последовательности векторов, сходящиеся к собственным векторам. Чаще всего последовательности собственных значений представляются как последовательности подобных матриц, сходящихся к треугольной или диагональной форме, что позволяет легко считывать собственные значения. Последовательности собственных векторов представляются соответствующими матрицами подобия. Метод Применяется к Производит Стоимость за шаг Сходимость Описание Алгоритм Ланчоса Эрмитовы m наибольшие/наименьшие собственные пары Итерация степеней Общая собственная пара с наибольшим значением O(n²) линейная Многократно применяет матрицу к произвольному начальному вектору и нормирует. Обратная итерация Общая собственная пара со значением, наиболее близким к μ линейная Итерация степеней для (A − μI)⁻¹ Итерация частного Рэлея Эрмитовы любая собственная пара кубическая Итерация степеней для (A − μiI)⁻¹, где μi для каждой итерации является частным Рэлея предыдущей итерации. Предварительная обратная итерация или алгоритм LOBPCG положительно определенные вещественно-симметричные собственные пары со значением, наиболее близким к μ. Обратная итерация с использованием предварительного решателя (приблизительный обратный к A). Метод бисекции Вещественно-симметричные тридиагональные любые линейные Использует метод бисекции для нахождения корней характеристического многочлена, поддерживаемый последовательностью Штурма. Итерация Лагерра Вещественно-симметричные тридиагональные любые кубические Использует метод Лагерра для нахождения корней характеристического многочлена, поддерживаемый последовательностью Штурма. QR-алгоритм Гессенова все собственные значения O(n²) кубический Разлагает A = QR, где Q ортогональна, а R треугольна, затем применяет следующую итерацию к RQ. все собственные пары 6n³ + O(n²) Алгоритм собственных значений Якоби Вещественно-симметричные все собственные значения O(n³) квадратичная Использует вращения Гивенса для попытки обнуления всех внедиагональных элементов. Это не удается, но усиливает диагональ. Метод "разделяй и властвуй" Эрмитовы тридиагональные все собственные значения O(n²) Разделяет матрицу на подматрицы, которые затем диагонализуются и рекомбинируются. все собственные пары Метод гомотопии Вещественно-симметричные тридиагональные все собственные пары O(n²) Строит вычислимый гомотопический путь из диагональной задачи собственных значений. Метод свернутого спектра Вещественно-симметричные собственные пары со значением, наиболее близким к μ. Предварительная обратная итерация, применяемая к (A − μI)² Алгоритм MRRR Вещественно-симметричные тридиагональные некоторые или все собственные пары O(n²) "Множественные относительно устойчивые представления" – выполняет обратную итерацию на LDLᵀ-разложении сдвинутой матрицы. Итерация Грама Общая собственная пара с наибольшим собственным значением сверхлинейная Многократно вычисляет грам-определитель и перемасштабирует детерминированно.

Iterative algorithms solve the eigenvalue problem by producing sequences that converge to the eigenvalues. Some algorithms also produce sequences of vectors that converge to the eigenvectors. Most commonly, the eigenvalue sequences are expressed as sequences of similar matrices which converge to a triangular or diagonal form, allowing the eigenvalues to be read easily. The eigenvector sequences are expressed as the corresponding similarity matrices. Method Applies to Produces Cost per step Convergence Description Lanczos algorithm Hermitian m largest/smallest eigenpairs Power iteration general eigenpair with largest value O(n^(2)) linear Repeatedly applies the matrix to an arbitrary starting vector and renormalizes. Inverse iteration general eigenpair with value closest to μ linear Power iteration for (A − μI)^(−1) Rayleigh quotient iteration Hermitian any eigenpair cubic Power iteration for (A − μiI)^(−1), where μi for each iteration is the Rayleigh quotient of the previous iteration. Preconditioned inverse iteration or LOBPCG algorithm positive definite real symmetric eigenpair with value closest to μ Inverse iteration using a preconditioner (an approximate inverse to A). Bisection method real symmetric tridiagonal any eigenvalue linear Uses the bisection method to find roots of the characteristic polynomial, supported by the Sturm sequence. Laguerre iteration real symmetric tridiagonal any eigenvalue cubic Uses Laguerre's method to find roots of the characteristic polynomial, supported by the Sturm sequence. QR algorithm Hessenberg all eigenvalues O(n^(2)) cubic Factors A = QR, where Q is orthogonal and R is triangular, then applies the next iteration to RQ. all eigenpairs 6n^(3) + O(n^(2)) Jacobi eigenvalue algorithm real symmetric all eigenvalues O(n^(3)) quadratic Uses Givens rotations to attempt clearing all off diagonal entries. This fails, but strengthens the diagonal. Divide and conquer Hermitian tridiagonal all eigenvalues O(n^(2)) Divides the matrix into submatrices that are diagonalized then recombined. all eigenpairs Homotopy method real symmetric tridiagonal all eigenpairs O(n^(2)) Constructs a computable homotopy path from a diagonal eigenvalue problem. Folded spectrum method real symmetric eigenpair with value closest to μ Preconditioned inverse iteration applied to (A − μI)^(2) MRRR algorithm real symmetric tridiagonal some or all eigenpairs O(n^(2)) "Multiple relatively robust representations" – performs inverse iteration on a LDLT decomposition of the shifted matrix. Gram iteration general Eigenpair with largest eigenvalue super linear Repeatedly computes the Gram product and rescales, deterministically.

Прямые расчеты

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

Треугольные матрицы

Поскольку определитель треугольной матрицы равен произведению её диагональных элементов, то если T — треугольная матрица, её собственные значения равны её диагональным элементам.

Эгенвекторы нормальных матриц 3х3

Если матрица 3×3 нормальна, то для нахождения собственных векторов можно использовать векторное произведение. Если λ является собственным значением матрицы A, то нулевое пространство матрицы (A - λI) перпендикулярно её столбцовому пространству. Векторное произведение двух независимых столбцов матрицы A будет лежать в нулевом пространстве. То есть, это будет собственный вектор, соответствующий λ. Поскольку столбцовое пространство в этом случае двумерно, собственное пространство должно быть одномерным, поэтому любой другой собственный вектор будет параллелен ему. Если матрица A не содержит двух независимых столбцов, но не является нулевой матрицей, векторное произведение всё равно можно использовать. В этом случае λ является собственным значением кратности 2, поэтому любой вектор, перпендикулярный столбцовому пространству, будет собственным вектором. Предположим, что v – ненулевой столбец матрицы A. Выберем произвольный вектор w, не параллельный v. Тогда (A - λI)v и (A - λI)w будут перпендикулярны v и, следовательно, будут собственными векторами матрицы A. Это не работает, когда матрица A не нормальна, так как нулевое пространство и столбцовое пространство не обязаны быть перпендикулярными для таких матриц.