Приближенный спектр большой матрицы

Я хочу вычислить спектр ( все собственные значения) большой разреженной матрицы (сотни тысяч строк). Это трудно.

Я готов согласиться на приближение. Существуют ли методы приближения для этого?

Хотя я надеюсь получить общий ответ на этот вопрос, я также был бы удовлетворен ответом в следующем конкретном случае. Моя матрица - нормализованный лапласиан большого графа. Собственные значения будут в диапазоне от 0 до 2 с большим количеством их сгруппированы около 1.

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

Чтобы получить приблизительную спектральную плотность, возьмите спектр ведущего подпространства Крылова размера 100 или около того и замените его дискретную плотность сглаженным вариантом.

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

Ответ Арнольда Ноймайера обсуждается более подробно в разделе 3.2 статьи «Приближение спектральных плотностей больших матриц» Лин Лин, Юсеф Саад и Чао Янг ​​(2016) .

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

Если вы не против думать о вещах, которые не являются собственными значениями, а функциями, которые в некотором смысле все еще говорят вам что-то о спектре, то я думаю, что вам стоит проверить некоторые работы Марка Эмбри в Университете Райса.

Вот еще один способ охарактеризовать спектр.

$\mathbf{A} v_k = \lambda_k v_k$$\mathbf{A}$

$$S(\omega) = \sum_k \frac{\pi^{-1}\sigma}{\sigma^2 + (\lambda_k - \omega)^2} = \frac{\sigma}{\pi} \mathop{\mathrm{Tr}} [\sigma^2 + (\omega - \mathbf{A})^2]^{-1}$$ $$S(\omega) = \frac{\sigma}{\pi}\langle z^T [\sigma^2 + (\omega - \mathbf{A})^2]^{-1} z \rangle$$ $z$$+1$$-1$$\sigma$$\omega$$[\sigma^2 + (\omega - \mathbf{A})^2]^{-1} z$$[\omega + i \sigma - \mathbf{A}]^{-1} [\omega - i \sigma - \mathbf{A}]^{-1}$чтобы минимизировать заполнение. Это позволяет оценить$S(\omega)$также для больших матриц. На практике кажется, что CG-решение не обязательно должно быть очень точным, и при вычислении среднего не требуется много векторов. Это может зависеть от проблемы.

Вышеуказанное, по-видимому, взвешивает части спектра более равномерно, чем аналогично смазанная спектральная плотность Крылова - попробуйте diag (linspace (0, 1, 150000)) - хотя, возможно, есть способ исправить это? Это несколько похоже на псевдоспектральный подход, но результат указывает (размытое) число собственных значений в окрестности точки$\omega$, а не обратное расстояние до ближайшего собственного значения.

РЕДАКТИРОВАТЬ : Лучшая альтернатива для вычисления вышеупомянутой величины состоит в том, чтобы вычислить чебышевские моменты (с помощью стохастической оценки, подобной описанной выше), а затем восстановить спектральную плотность из них. Это не требует ни матричных инверсий, ни отдельных вычислений для каждого$\omega$. See http://theorie2.physik.uni-greifswald.de/downloads/publications/LNP_chapter19.pdf and references therein.

See the paper "On Sampling-based Approximate Spectral Decomposition" by Sanjiv Kumar, Mehryar Mohri & Ameet Talwalkar (ICML 2009.). It uses sampling of columns of your matrix.

Since your matrix is symmetric you should do the following:

Let A be your n*n matrix. You want to reduce the computation of the eigenvalues of an n*n matrix to the computation of the eigenvalues of an k*k matrix. First choose your value of k. Let's say you choose k=500, since you can easily compute the eigenvalues of a 500*500 matrix. Then, randomly choose k columns of the matrix A. Contruct the matrix B that keeps only these columns, and the corresponding rows.

B = A(x,x) for a random set of k indexes x

B is now a k*k matrix. Compute the eigenvalues of B, and multiply them by (n/k). You now have k values which are approximately distributed like the n eigenvalues of A. Note that you get only k values, not n, but their distribution will be correct (up to the fact that they are an approximation).

You can always use the Gershgorin circle Theorem bounds to approximate the eigenvalues.

If the off-diagonal terms are small, the diagonal itself is a good approximation of the spectrum. Otherwise if you end up with an approximation of the eigenspace (by other methods) you could try to express the diagonal entries in this system. This will lead to a matrix with smaller off-diagonal terms and the new diagonal will be a better approximation of the spectrum.