Ctrl+K

Матричные разложения

Матричные разложения#

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

Сингулярное разложение (SVD)#

original image

  • \(r=\operatorname{rank}(A)\)

  • \(U, V\) - унитарные (\(U^{*}U=E\))

  • \(\sigma_1 \geq \ldots \geq \sigma_r>0\) - ненулевые сингулярные числа. Для эрмитовых матриц равны модулям собственных чисел.

  • Столбцы \(U, V\) - сингулярные вектора.

  • Существует всегда.

Note: SVD также может быть определена для \(U \in \mathbb{C}^{m \times p}\), \(\Sigma \in \mathbb{R}^{p \times p}\) и \(V \in \mathbb{C}^{n \times p}, p=\min \{n . m\}\)

Алгоритмы построения#

  • Спектральное разложение \(A A^*\) или \(A^* A\) совпадает со спектральным разложением (собственное разложение).

  • Стабильный алгоритм \(\mathcal{O}\left(m n^2\right)\) , \((m>n):\)

  1. Bidiagonalize \(A\) by Householder reflections

\[ A=U_1 B V_1^* \]

  1. Find SVD of \(B=U_2 \Sigma V_2^*\) by spectral decomposition of \(T\) ( 2 options).

  2. \(U=U_1 U_2, \quad V=V_1 V_2\)

Области применения#

  • Сжатие данных

  • Эффективное вычисление псевдообратной матрицы

  • Определение ключевых признаков в Data Science

Спектральное разложение#

original image

  • \(\lambda_1, \ldots, \lambda_n\) - собственные числа.

  • Столбцы \(S\) - собственные вектора.

  • Не всегда существует (в отличие от верхнетреугольного вида)

Алгоритмы построения#

  • Если \(A=A^*\), Jacobi method: \(\mathcal{O}\left(n^3\right)\)

  • Если \(A A^*=A^* A\), QR algorithm: \(\mathcal{O}\left(n^3\right)\)

  • Если \(A A^* \neq A^* A\) ассимптотика также \(\mathcal{O}\left(n^3\right)\)

Области применения#

  • Полное спектральное разложение используется редко, т.к. для вычисления необходимы все собственные вектора.

  • Для вычисления собственных чисел лучше использовать верхнетреугольную форму (существует алгоритм за \(\mathcal{O}\left(n^3\right)\))

LU-разложение (Разложение Холецкого)#

original image

  • Существует только при ненулевых лидирующих минорах.

  • Уникально.

Если матрица \(A\) - симметрична и положительно определена, LU-разложение можно модифицировать в разложение Холецкого: original image

Алгоритмы построения#

  • Different versions of Gaussian elimination, \(\mathcal{O}\left(n^3\right)\) flops. In LU for stability use permutation of rows or columns (LUP).

  • \(\mathcal{O}\left(n^3\right)\) can be decreased for sparse matrices by appropriate permutations, e.g.

    • minimum degree ordering

    • Cuthill-Mckee algorithm

  • Banded matrix with bandwidth \(b\)

original image

can be decomposed using \(\mathcal{O}\left(n b^2\right)\) flops.

Области применения#

LU и более высокоуровневые модификации LDL, Холецкий используются для

  • Решения СЛАУ. Пусть \(A=L U\), тогда сложность решения \(A x=b\) есть \(\mathcal{O}\left(n^2\right)\) :

    1. Решим простую систему: \(L y=b\)

    2. Решим обратную систему: \(U x=y\)

  • Обращения матриц

  • Вычисления детерминанта

  • Как часть других разложений.

Пример 1 - решить СЛАУ с помощью LU-разложения#

Пусть СЛАУ:

\[ Ax=b \]

Тогда хотим:

\[ Ax = LUx = b \]

Для чего последовательно решаем:

\[ Ly=b \]
\[ Ux=y \]

Будем использовать следующие обозначения для элементов матриц: \(A=\left(a_{i j}\right), L=\left(l_{i j}\right), U=\left(u_{i j}\right), i, j=1 \ldots n\); причём диагональные элементы матрицы \(L: l_{i i}=1, i=1 \ldots n\).

Найти матрицы \(L\) и \(U\) можно следующим образом (выполнять шаги следует строго по порядку, так как следующие элементы находятся с использованием предыдущих):

  1. Цикл і от 1 до n

    1. Цикл ј от 1 до \(\mathrm{n}\)

      1. \(u_{i j}=0, l_{i j}=0\)

      2. \(l_{i i}=1\)

  2. Цикл і от 1 до \(\mathrm{n}\)

    1. Цикл ј от 1 до \(\mathrm{n}\)

      1. Если \(i \leq j: u_{i j}=a_{i j}-\sum_{k=1}^{i-1} l_{i k} \cdot u_{k j}\)

      2. Если \(i>j: l_{i j}=\left(a_{i j}-\sum_{k=1}^{j-1} l_{i k} \cdot u_{k j}\right) / u_{j j}\)

В итоге мы получим матрицы - \(L\) и \(U\).

import numpy as np

def lu_decomposition(A):
    n = A.shape[0]
    U = np.zeros((n, n))
    L = np.eye(n)

    for i in range(n):
        U[i, i:] = A[i, i:] - L[i, :i] @ U[:i, i:]
        L[i+1:, i] = (A[i+1:, i] - L[i+1:, :i] @ U[:i, i]) / U[i, i]

    return L, U

A = np.array([[1, 7, 3], [4, 5, 6], [7, 8, 9.]])

L, U = lu_decomposition(A)

A - L @ U

array([[ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
       [ 0.00000000e+00, -1.77635684e-15,  0.00000000e+00]])

b = np.array([1., 2., 3.]) # столбец свободных членов

# Решение Ly = b

y = b*0

for i in range(L.shape[0]):
    y[i] = b[i] - L[i, :i] @ y[:i]

# Решение Ux = y

x = b*0

for i in range(U.shape[0]-1, -1, -1):
    x[i] = (y[i] - U[i, i+1:] @ x[i+1:])/U[i, i]

# Проверка

x - np.linalg.inv(A) @ b

array([ 1.11022302e-16, -8.32667268e-17,  1.66533454e-16])

Пример 2 - решение МНК через Холецкого#

Вспомним что для заданной матрицы объекты-признаки \(X\) наилучшие веса в MSE находились по формуле:

\[ \mathbf{w} = (X^\top X)^{-1} \cdot X^\top \cdot \mathbf{y} \]

Что эквивалентно решению системы:

\[ (X^\top X) \cdot \mathbf{w} = X^\top \cdot \mathbf{y} \]
\[ A \cdot \mathbf{w} = \mathbf{b} \]

С симметричной положительно определенной матрицей \(A=(X^\top X)\) и столбцом свободных членов \(\mathbf{b} = X^\top \cdot \mathbf{y}\). Намечается расзложение Холецкого, для тех кто не понял :).

Его плюс в данном случае - скорость выполнения (т.к. нам не нужно явно считать обратную).

import numpy as np
from numba import njit

def cholesky_solve(A, b):
    L = np.linalg.cholesky(A)

    # Решение Ly = b
    y = np.zeros_like(b)
    for i in range(L.shape[0]):
        y[i] = (b[i] - L[i, :i] @ y[:i])/L[i, i]

    # Решение L.T x = y
    x = np.zeros_like(b)
    for i in range(L.shape[1]-1, -1, -1):
        x[i] = (y[i] - L[i+1:, i] @ x[i+1:])/L[i, i]

    return x

A = np.array([[1, 7, 3], [4, 5, 6], [7, 8, 9.]])

try:
    L = np.linalg.cholesky(A) # Только для положительно определенных!
except Exception as error:
    print(error)

Matrix is not positive definite

# Фит многочленами с прошлого семинара

import numpy as np

# Генерируем выборку для исследования

x = 2*np.random.random(1000000).reshape((-1, 1)) - 1
noise = 0.1*np.random.normal(size = x.shape)
y_label = -4 + 3*x - 2*x**2 + 5*x**3 + noise

# Напишем функцию для создания матрицы объекты-признаки
def create_X(x, m):
    X = np.ones(x.shape) # Сразу же добавим intercept
    for i in range(1, m + 1):
        X = np.concatenate([X, x**i], axis = 1)
    return X

X = create_X(x, 5)
A = X.T @ X

L = np.linalg.cholesky(A)

L @ L.T - A

array([[ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
       [ 0.00000000e+00,  5.82076609e-11,  0.00000000e+00,
        -2.91038305e-11,  0.00000000e+00, -2.91038305e-11],
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
       [ 0.00000000e+00, -2.91038305e-11,  0.00000000e+00,
         0.00000000e+00,  2.84217094e-14,  0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         2.84217094e-14,  0.00000000e+00,  0.00000000e+00],
       [ 0.00000000e+00, -2.91038305e-11,  0.00000000e+00,
         0.00000000e+00,  0.00000000e+00,  0.00000000e+00]])

%%time
w_pinv = np.linalg.inv(A) @ X.T @ y_label # Вычисление обратной O(n^3)

CPU times: user 0 ns, sys: 15.6 ms, total: 15.6 ms
Wall time: 55.1 ms

%%time
w_chol = cholesky_solve(A, X.T @ y_label) # А решение СЛАУ за O(n^2)

CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 5.23 ms

w_chol - w_pinv

array([[ 7.63833441e-14],
       [-2.41584530e-13],
       [-6.11954931e-13],
       [ 6.39488462e-13],
       [ 7.05184609e-13],
       [-5.01479500e-13]])
Содержание