Собственные векторы треугольных матриц: быстрый путь через numpy.linalg.eig

Ускорьте вычисление собственных векторов для верхнетреугольных матриц: используйте numpy.linalg.eig. Бенчмарки, код и разбор, почему это быстрее и надёжнее.

Когда на входе треугольная матрица, собственные значения «считываются» прямо с диагонали. Возникает желание обойти универсальные процедуры и на скорую руку написать что-то быстрее для диагонализации или извлечения собственных векторов. Подвох в том, что не всякая треугольная матрица диагонализуема, а если нужны только собственные векторы, готовый путь уже может быть самым быстрым.

Постановка задачи

Требуется вычислять собственные векторы для множества треугольных матриц без лишних накладных расходов. Диагональ сразу дает собственные значения, поэтому решение (A − λI)x = 0 для каждого λ выглядит заманчиво. Но если у вас сотни матриц и примерно столько же собственных значений на каждую, отдельные решения для каждого λ быстро накапливают время. Практический вопрос в том, не «слишком ли общий» здесь numpy.linalg.eig или он уже оптимизирован под треугольный ввод.

Пример кода: сравнение подходов по скорости

Ниже — фрагмент, который сопоставляет производительность собственных разложений с QR- и LU-факторизациями на плотных случайных матрицах и их верхнетреугольных версиях.

import numpy as np
import scipy.linalg as sla
np.random.seed(42)

for dim in [100, 1000]:
    print(f"benchmark on size {dim}")
    mat = np.random.normal(size=(dim, dim))
    print("dense input")
    %timeit np.linalg.eig(mat)
    %timeit np.linalg.qr(mat)
    %timeit sla.lu_factor(mat)
    print("upper triangular input")
    tri = np.random.normal(size=(dim, dim))
    tri = np.triu(tri)
    %timeit np.linalg.eig(tri)
    %timeit np.linalg.qr(tri)
    %timeit sla.lu_factor(tri)

Что происходит на самом деле

Эффективные алгоритмы для собственных векторов обычно одновременно получают и собственные значения, поэтому, вызывая единую процедуру, вы ничего не теряете. На практике numpy.linalg.eig хорошо настроен даже для уже верхнетреугольных входов. Изучение реализации показывает, что OpenBLAS выполняет QR-факторизацию, чтобы перейти к верхнетреугольной форме, а затем использует рутину LAPACK strevc3_ для получения собственных значений и векторов. Если на вход сразу подать верхнетреугольную матрицу, шаг QR особенно быстрый.

Эмпирические замеры это подтверждают. На верхнетреугольном вводе размера 100 np.linalg.eig укладывается примерно в сотни микросекунд, а для размера 1000 — в сотни миллисекунд. В тех же условиях запуск eig на полной плотной матрице занимает около 11.7 ms для размера 100 и примерно 1.31 s для 1000. Иными словами, на треугольных входах достигается заметное ускорение по сравнению с плотными.

Решение: просто используйте numpy.linalg.eig для треугольной матрицы

Если у вас уже есть верхнетреугольная матрица, вызывайте eig напрямую. Вы избегаете накладных расходов на решения для каждого λ и используете быстрый путь. Минимальный пример:

import numpy as np
np.random.seed(42)

n = 1000
tri = np.triu(np.random.normal(size=(n, n)))
vals, vecs = np.linalg.eig(tri)

В более широких испытаниях этот подход оказывался быстрее примерно в 10–20 раз для верхнетреугольных матриц по сравнению с обработкой плотных матриц того же размера. Ранний бенчмарк показывает ту же картину: eig на треугольных входах завершается примерно за 438 µs против 11.7 ms для размера 100 и около 157 ms против 1.31 s для 1000.

Если вы подумываете напрямую вызвать strevc3_ из Python, простого пути нет. Хотя внутри SciPy и NumPy присутствуют процедуры strevc3_ или dtrevc3, ни одна библиотека не предоставляет к ним открытой Python-обертки.

Важное замечание о диагонализации

Существенно отличать вычисление собственных векторов от диагонализации. Не каждая треугольная матрица диагонализуема. Простой контрпример — [[1, 1], [0, 1]]: у нее одно собственное значение и недостаточно линейно независимых собственных векторов, чтобы образовать базис, поэтому диагонализовать ее нельзя. Если ваша конечная цель — именно диагонализация, убедитесь, что она теоретически возможна; иначе попытка собрать диагональный вид провалится, как бы быстро вы ни считали собственные векторы.

Почему это важно

При обработке множества матриц сокращение постоянных факторов радикально влияет на время выполнения. Оптимизированное собственное разложение использует быстрые пути для треугольных входов и избавляет от повторных решений для каждого λ, которые добавляют накладные расходы. Код остаётся короче и проще в сопровождении — без потери производительности.

Выводы

Если матрица верхнетреугольная, вызывайте numpy.linalg.eig и передавайте данные как есть. Не тратьте ресурсы на отдельные решения для каждого собственного значения. Помните: «треугольная» не означает «диагонализуемая»; пример [[1, 1], [0, 1]] показывает, что диагонализация может быть невозможна, даже когда собственные значения легко читаются с диагонали. Если вы рассчитывали подключить strevc3_ из Python, учтите: хотя рутина есть внутри, в NumPy и SciPy нет простого публичного обертчика.