Самый быстрый способ решения линейных наименьших квадратов

Question

В https://math.stackexchange.com/a/2233298/340174 упоминается, что решение линейного уравнения "M · x = b" (матрица M является квадратной) является медленным, если выполняется посредством разложения LU (но еще медленнее, используя разложение QR).Теперь я заметил, что numpy.linalg.solve фактически использует разложение LU.По правде говоря, я хочу решить «V · x = b» для не квадратной матрицы Вандермонда V для наименьших квадратов.Я не хочу регуляризации.Я вижу несколько подходов:

1. Решите "V · x = b" с numpy.linalg.lstsq, который использует Fortran "xGELSD" на основе SVD.SVD должен быть даже медленнее, чем разложение LU, но мне не нужно вычислять «(V ^ T · V)».
2. Solve »(V ^ T · V) · x = (V ^ T· Б) "с numpy.linalg.solve, который использует разложение LU.
3. Решите" A · x = b "с numpy.linalg.solve, который использует разложение LU, но вычисляет" A = xV ^ T · V "напрямуюв соответствии с https://math.stackexchange.com/a/3155891/340174

В качестве альтернативы я мог бы использовать новейший solve от scipy (https://docs.scipy.org/doc/scipy-1.2.1/reference/generated/scipy.linalg.solve.html), который может использовать диагональное вращение для симметричной матрицы "A" (что быстрее, чем использованиеРазложение LU, я думаю), но мой scipy застрял на 1.1.0, поэтому у меня нет доступа к этому.

С https://stackoverflow.com/a/45535523/4533188 кажется, что solve быстрее, чем lstsq, включая вычисление «V ^ T · V», но когда я попробовал, lstsq был быстрее. Может быть, я делаю что-то не так?

Какой самый быстрый способ решения моей линейнойпроблема?

---

Нет реальных опций

• statsmodels.regression.linear_model.OLS.fit используется либо псевдообратная Мура-Пенроуза, либо QR-факторизация + np.linalg.inv + np.linalg.svd + numpy.linalg.solve, что немне не кажется слишком эффективным.
• sklearn.linear_model.LinearRegression использует scipy.linalg.lstsq.
• scipy.linalg.lstsq использует также xGELSD.
• Я ожидаю вычисления обратного значения "(V ^ T · V) "быть довольно дорогим, поэтому я отбросил прямое вычисление" x = (V ^ T · V) ^ - 1 · (V ^ T · b) "

Answer 1

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

Я думаю, что некоторые вещи могут быть здесь перепутаны, поэтому я выделю следующее:

1. Точное решение V x = b с использованием LU
2. Точное решение V x = b с использованием QR
3. Решение наименьших квадратов V x = b с использованием QR
4. Решение наименьших квадратов V x = b с использованием SVD
5. Точное решение V^T V x = V^T b с использованием LU
6. Точное решение V^T V x = V^T b с использованием Cholesky

Первый ответ maths.stackexchange, на который вы ссылались, относится к случаям 1 и 2. Когда он говорит, что LU медленный, это означает относительно методов для конкретных типов матриц, например, положительно определенный, треугольный, полосатый, ...

Но я думаю, что вы на самом деле спрашиваете о 3-6. Последняя ссылка на стек-поток утверждает, что 6 быстрее, чем 4. Как вы сказали, 4 должно быть медленнее, чем 3, но 4 - единственная, которая работает для ранг-дефицитных V. 6 должно быть быстрее, чем 5 в целом.

Мы должны убедиться, что вы сделали 6, а не 5. Чтобы использовать 6, вам нужно будет использовать scipy.linalg.solve с assume_a="pos". В противном случае вы бы сделали 5.

Я не нашел ни одной функции, которая выполняет 3 в numpy или scipy. Подпрограмма Lapack - это xGELS, которая, похоже, не раскрывается в scipy. Вы должны быть в состоянии сделать это с помощью scupy.linalg.qr_multiply, за которым следует scipy.linalg.solve_triangular.

Answer 2

Попробуйте scipy.linalg.lstsq(), используя lapack_driver='gelsy'!

Давайте рассмотрим различные процедуры для решения линейного метода наименьших квадратов и приближений:

• numpy.linalg.lstsq() упаковывает LAPACK xGELSD(), как показано в umath_linalg.c.src в строке 2841+.Эта подпрограмма переводит матрицу V в бидиагональную форму, используя стратегию разделения и завоевания, и вычисляет SVD этой бидиагональной матрицы.

• scipy's scipy.linalg.lstsq() обертывает LAPACK xGELSD(), xGELSY() и xGELSS(): аргумент lapack_driver может быть изменен для переключения с одного на другой.Согласно тесту LAPACK, xGELSS() медленнее, чем xGELSD(), а xGELSY() примерно на 5 быстрее, чем xGELSD().xGELSY() использует QR-факторизацию V с поворотом столбца.И хорошая новость в том, что этот переключатель был уже доступен в Scipy 1.1.0 !

• В LAPACK xGELS() используется QR-разложение матрицы V, но предполагается, что эта матрица имеет полный ранг.Согласно тесту LAPACK, можно ожидать, что dgels() будет примерно в 5 раз быстрее, чем dgelsd(), но он также более уязвим для номера условия матрицы и может стать неточным.Подробности и дальнейшие ссылки см. В Разница между результатами C ++ (LAPACK, sgels) и Python (Numpy, lstsq) .xGELS () доступен в cython-lapack интерфейсе scipy .

Хотя очень заманчиво, вычисление и использование V^T·V для решения нормального уравнения, скорее всего, не лучший путь.Действительно, точность подвергается опасности из-за номера условия этой матрицы, около квадрат номера условия матрицы V .Поскольку матрицы Вандермонда имеют тенденцию быть плохо обусловленными, за исключением матриц дискретного преобразования Фурье , это может стать опасным ... Наконец, вы даже можете продолжать использовать xGELSD(), чтобы избежать проблемсвязанные с кондиционированием.Если вы переключитесь на xGELSY(), рассмотрите , оценивая ошибку .