SVD e aproximações de matrizes

1.1–1.2 - Brunton: Data-Driven Science and Engineering

Prof. Dr. Raphael Teixeira

Universidade Federal do Pará

Muitos dados - Poucos padrões

Sistemas Complexos:

• Muitos dados: séries temporais, imagens;
• Poucos padrões (posto) dominantes: modos, polos, frequências;

Ex: identificação de sistemas

Centenas de amostras, mas apenas dois padrões dominantes: os padrões são os seus dois polos.

A matriz de dados

Cada coluna de \(\mathbf{X} = [\mathbf{x}_1 \ \cdots \ \mathbf{x}_m] \in \mathbb{C}^{n \times m}\) é um conjunto de medições: uma imagem, ou amostras de entrada e saída de um sistema.

• As colunas são snapshots \(\mathbf{x}_k = \mathbf{x}(k\Delta t)\); em geral \(n \gg m\): uma matriz alta e estreita.

Uma matriz como soma de padrões

Cada padrão é uma coluna \(\mathbf{u}_k\) vezes uma linha \(\mathbf{v}_k^T\), com um peso \(\sigma_k\) (valor singular).

Da soma de padrões à SVD

Para \(\mathbf{X} \in \mathbb{C}^{n \times m}\) com \(n \ge m\), a matriz é uma soma de até \(m\) padrões:

\[ \mathbf{X} = \sigma_1 \mathbf{u}_1 \mathbf{v}_1^* + \sigma_2 \mathbf{u}_2 \mathbf{v}_2^* + \cdots + \sigma_m \mathbf{u}_m \mathbf{v}_m^* = \sum_{k=1}^{m} \sigma_k \mathbf{u}_k \mathbf{v}_k^* \qquad (1.4) \]

\[ \color{#1B365D}{\Large \mathbf{X} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^*} \]

\[ \color{#B23A48}{\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_m \ge 0} \qquad \text{os padrões vêm ranqueados por importância} \]

A SVD completa: \(\mathbf{U}\), \(\boldsymbol{\Sigma}\) e \(\mathbf{V}\)

\[ \boldsymbol{\Sigma} = \left[\begin{array}{ccc} \sigma_1 & & \\ & \ddots & \\ & & \sigma_m \\ \hline & \mathbf{0} & \end{array}\right] \begin{array}{l} \left.\vphantom{\begin{array}{c} \sigma_1 \\ \ddots \\ \sigma_m \end{array}}\right\} m \times m \ \text{(diagonal)} \\ \left.\vphantom{\begin{array}{c} \mathbf{0} \end{array}}\right\} (n-m) \times m \end{array} \]

• \(\mathbf{U}\) (\(n \times n\)): suas colunas \(\mathbf{u}_k\) são os padrões ao longo das colunas de \(\mathbf{X}\).
• \(\boldsymbol{\Sigma}\) (\(n \times m\)): diagonal com os valores singulares de \(\mathbf{X}\), \(\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_m \ge 0\) — os pesos dos padrões;
• \(\mathbf{V}^*\) (\(m \times m\)): suas linhas \(\mathbf{v}_k^*\) são os padrões ao longo das linhas de \(\mathbf{X}\).

As colunas de \(\mathbf{U}\) além da \(m\)-ésima (\(\mathbf{\hat{U}}^\perp\)) multiplicam só zeros de \(\boldsymbol{\Sigma}\).

Diferente da FFT, que usa uma base fixa (senoides), as bases \(\mathbf{U}\) e \(\mathbf{V}\) da SVD são extraídas dos próprios dados, sem conhecimento prévio do sistema.

A decomposição SVD

Definição

Toda matriz \(\mathbf{X} \in \mathbb{C}^{n \times m}\) pode ser escrita como \[ \mathbf{X} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^* \] com \(\mathbf{U}\) (\(n \times n\)) e \(\mathbf{V}\) (\(m \times m\)) unitárias e \(\boldsymbol{\Sigma}\) (\(n \times m\)) diagonal, \(\sigma_k \ge 0\).

• Posto de \(\mathbf{X}\) = número de \(\sigma_k \neq 0\).
• Se \(n \ge m\), só as \(m\) primeiras colunas de \(\mathbf{U}\) multiplicam valores não nulos: \[ \mathbf{X} = \mathbf{\hat{U}}\, \boldsymbol{\hat{\Sigma}}\, \mathbf{V}^* \] a SVD econômica.

Computando a SVD

Na prática, basta uma linha: a SVD vem pronta do LAPACK, base do NumPy e do MATLAB.

• Python
• MATLAB

import numpy as np

X = np.array([[3,2,2], [2,3,-2], [1,0,1], [0,1,1], [2,2,0]])   # dados 5 x 3
U, S, VT = np.linalg.svd(X, full_matrices=False)               # SVD econômica

X = [3 2 2; 2 3 -2; 1 0 1; 0 1 1; 2 2 0];   % dados 5 x 3

[U,S,V] = svd(X,'econ');                     % SVD econômica

\[ \overset{\displaystyle \mathbf{X}\ (5 \times 3)}{\begin{bmatrix} 3 & 2 & 2 \\ 2 & 3 & -2 \\ 1 & 0 & 1 \\ 0 & 1 & 1 \\ 2 & 2 & 0 \end{bmatrix}} = \overset{\displaystyle \mathbf{\hat{U}}\ (5 \times 3)}{\begin{bmatrix} -0.63 & -0.58 & -0.11 \\ -0.58 & \phantom{-}0.71 & \phantom{-}0.02 \\ -0.13 & -0.34 & -0.36 \\ -0.13 & -0.21 & \phantom{-}0.93 \\ -0.48 & \phantom{-}0.05 & -0.04 \end{bmatrix}} \raise{1.25em}{\overset{\displaystyle \boldsymbol{\hat{\Sigma}}\ (3 \times 3)}{\begin{bmatrix} \color{#B23A48}{5.84} & 0 & 0 \\ 0 & \color{#B23A48}{3.29} & 0 \\ 0 & 0 & \color{#B23A48}{1.05} \end{bmatrix}}} \raise{1.25em}{\overset{\displaystyle \mathbf{V}^T\ (3 \times 3)}{\begin{bmatrix} -0.71 & -0.70 & -0.06 \\ -0.17 & \phantom{-}0.26 & -0.95 \\ -0.68 & \phantom{-}0.66 & \phantom{-}0.30 \end{bmatrix}}} \]

O NumPy devolve os valores singulares como vetor e \(\mathbf{V}^T\); o MATLAB devolve \(\boldsymbol{\Sigma}\) como matriz e \(\mathbf{V}\).

Soma diádica

Como \(\boldsymbol{\Sigma}\) é diagonal, \(\mathbf{X}\) é uma soma de matrizes de posto 1, em hierarquia:

\[ \mathbf{X} = \sum_{k=1}^{m} \sigma_k \mathbf{u}_k \mathbf{v}_k^* = \sigma_1 \mathbf{u}_1 \mathbf{v}_1^* + \sigma_2 \mathbf{u}_2 \mathbf{v}_2^* + \cdots + \sigma_m \mathbf{u}_m \mathbf{v}_m^* \]

com \(\mathbf{u}_k\) e \(\mathbf{v}_k\) as \(k\)-ésimas colunas de \(\mathbf{U}\) e \(\mathbf{V}\).

• \(\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_m \ge 0\): termos ranqueados, cada um menos importante que o anterior.
• Em muitos sistemas os \(\sigma_k\) decaem rapidamente, e basta truncar a soma em um posto \(r\).

SVD truncada: aproximação matricial

\[ \color{#1B365D}{\mathbf{X} \approx \mathbf{\tilde{X}} = \sum_{k=1}^{r} \sigma_k \mathbf{u}_k \mathbf{v}_k^* = \mathbf{\tilde{U}} \boldsymbol{\tilde{\Sigma}} \mathbf{\tilde{V}}^*} \]

• \(\mathbf{\tilde{U}}\), \(\mathbf{\tilde{V}}\): as \(r\) primeiras colunas de \(\mathbf{U}\), \(\mathbf{V}\); \(\boldsymbol{\tilde{\Sigma}}\): o bloco \(r \times r\).
• Se \(r <\) posto de \(\mathbf{X}\), \(\mathbf{\tilde{X}}\) só aproxima \(\mathbf{X}\).
• \(\mathbf{\tilde{U}}\) leva das medições (dimensão \(n\)) aos padrões (dimensão \(r\)).

Teorema de Eckart–Young

Resultado de Schmidt (o de Gram–Schmidt), redescoberto por Eckart e Young:

Teorema 1.1 (Eckart–Young)

A melhor aproximação de posto \(r\) de \(\mathbf{X}\), no sentido de mínimos quadrados, é a SVD truncada de posto \(r\): \[ \underset{\mathbf{\tilde{X}},\ \text{s.a.}\ \operatorname{posto}(\mathbf{\tilde{X}}) = r}{\operatorname{argmin}} \ \big\| \mathbf{X} - \mathbf{\tilde{X}} \big\|_F = \mathbf{\tilde{U}} \boldsymbol{\tilde{\Sigma}} \mathbf{\tilde{V}}^* \]

ou seja, \[ \mathbf{X} \approx \mathbf{\tilde{X}} = \sum_{k=1}^{r} \sigma_k \mathbf{u}_k \mathbf{v}_k^* = \sigma_1 \mathbf{u}_1 \mathbf{v}_1^* + \sigma_2 \mathbf{u}_2 \mathbf{v}_2^* + \cdots + \sigma_r \mathbf{u}_r \mathbf{v}_r^* \]

Norma de Frobenius: \(\|\mathbf{X}\|_F = \big( \sum_{i=1}^{n} \sum_{j=1}^{m} |X_{ij}|^2 \big)^{1/2}\), a norma 2 da matriz “vetorizada” \(\mathbf{X}(:)\).

Erro da aproximação na norma de Frobenius

O erro da SVD truncada de posto \(r\) é conhecido exatamente, e nenhuma outra matriz de posto \(r\) erra menos:

\[ \big\| \mathbf{X} - \mathbf{\tilde{X}} \big\|_F^2 = \sum_{k=r+1}^{m} \sigma_k^2. \qquad (1.7) \]

Por quê?

O resíduo é a soma dos termos descartados, \(\mathbf{X} - \mathbf{\tilde{X}} = \sum_{k=r+1}^{m} \sigma_k \mathbf{u}_k \mathbf{v}_k^* = \mathbf{U}_{\text{rem}} \boldsymbol{\Sigma}_{\text{rem}} \mathbf{V}_{\text{rem}}^*\). Como \(\|\mathbf{A}\|_F^2 = \operatorname{tr}(\mathbf{A}^*\mathbf{A})\) e \(\mathbf{U}_{\text{rem}}^*\mathbf{U}_{\text{rem}} = \mathbf{V}_{\text{rem}}^*\mathbf{V}_{\text{rem}} = \mathbf{I}\), as bases ortonormais não alteram a norma: \(\|\mathbf{X} - \mathbf{\tilde{X}}\|_F^2 = \|\boldsymbol{\Sigma}_{\text{rem}}\|_F^2 = \sigma_{r+1}^2 + \cdots + \sigma_m^2\).

Como esse erro escala com a magnitude de \(\mathbf{X}\), costuma ser mais útil o erro relativo \(\;\big\| \mathbf{X} - \mathbf{\tilde{X}} \big\|_F^2 \,/\, \|\mathbf{X}\|_F^2 \;\) (1.8):

• Sinais medidos (tensões, correntes): (1.8) é a fração da energia que falta (\(\|\mathbf{X}\|_F^2\) = soma dos quadrados das amostras).
• Dados com média subtraída: é a fração da variância que falta (PCA, §1.5).

Aproximação ótima na norma 2

A SVD truncada também é ótima na norma 2 (norma espectral), induzida pela norma 2 de vetores:

\[ \underset{\mathbf{\tilde{X}},\ \text{s.a.}\ \operatorname{posto}(\mathbf{\tilde{X}}) = r}{\operatorname{argmin}} \ \big\| \mathbf{X} - \mathbf{\tilde{X}} \big\|_2 = \mathbf{\tilde{U}} \boldsymbol{\tilde{\Sigma}} \mathbf{\tilde{V}}^*, \qquad (1.9) \qquad \|\mathbf{X}\|_2 = \max_{\mathbf{v} \neq \mathbf{0}} \frac{\|\mathbf{X}\mathbf{v}\|_2}{\|\mathbf{v}\|_2}. \]

Nessa norma, o erro é ainda mais simples: \[ \big\| \mathbf{X} - \mathbf{\tilde{X}} \big\|_2 = \sigma_{r+1}. \qquad (1.10) \]

Basta expandir o resíduo, \[ \mathbf{X} - \mathbf{\tilde{X}} = \sum_{k=r+1}^{m} \sigma_k \mathbf{u}_k \mathbf{v}_k^*, \qquad (1.11) \] e usar a ortonormalidade dos \(\mathbf{v}_k\): o máximo de \(\|(\mathbf{X} - \mathbf{\tilde{X}})\mathbf{v}\|_2\) sobre vetores unitários é \(\sigma_{r+1}\), atingido em \(\mathbf{v} = \mathbf{v}_{r+1}\).

Verificação numérica das expressões de erro

• Python
• MATLAB

rng = np.random.default_rng(0)
X = rng.standard_normal((200, 50))
U, S, VT = np.linalg.svd(X, full_matrices=False)

r = 10
Xr = U[:, :r] @ np.diag(S[:r]) @ VT[:r, :]         # SVD truncada de posto r (1.5)

print(np.linalg.norm(X - Xr, 'fro')**2, np.sum(S[r:]**2))   # (1.7)
print(np.linalg.norm(X - Xr, 2), S[r])                      # (1.10)

6438.2718056637295 6438.2718056637295
17.205369539578513 17.205369539578506

rng(0); X = randn(200,50);
[U,S,V] = svd(X,'econ');
s = diag(S);                                % valores singulares (vetor)

r = 10;
Xr = U(:,1:r)*S(1:r,1:r)*V(:,1:r)';         % SVD truncada de posto r (1.5)

disp([norm(X - Xr,'fro')^2, sum(s(r+1:end).^2)])   % (1.7)
disp([norm(X - Xr, 2), s(r+1)])                    % (1.10)

Os índices em Python começam em 0: S[r] é \(\sigma_{r+1}\) e S[r:] são \(\sigma_{r+1}, \dots, \sigma_m\); no MATLAB, s(r+1) e s(r+1:end).

Exemplo: compressão de imagem

Uma imagem em tons de cinza é uma matriz real \(\mathbf{X} \in \mathbb{R}^{n \times m}\) (\(n\) pixels na vertical, \(m\) na horizontal). Aqui, a mesma foto do livro: o cachorro Mordecai (\(2000 \times 1500\) pixels).

• Python
• MATLAB

A = plt.imread("imgs/dog.jpg")                    # foto colorida (RGB)
X = np.mean(A, axis=-1)                           # RGB -> tons de cinza
n, m = X.shape
U, S, VT = np.linalg.svd(X, full_matrices=False)  # SVD econômica

aproxima = lambda r: U[:, :r] @ np.diag(S[:r]) @ VT[:r, :]  # SVD truncada (1.5)
armazenamento = lambda r: r * (n + m + 1) / (n * m)        # fração guardada

A = imread('imgs/dog.jpg');              % foto colorida (RGB)
X = mean(double(A), 3);                  % RGB -> tons de cinza
[n,m] = size(X);
[U,S,V] = svd(X,'econ');  s = diag(S);   % SVD econômica
for r = [5 20 100]                       % SVD truncada (1.5) para vários postos
    Xapprox = U(:,1:r)*S(1:r,1:r)*V(:,1:r)';
    figure, imagesc(Xapprox), axis off, colormap gray
    title(sprintf('r = %d  (%.1f%% dos dados)', r, 100*r*(n+m+1)/(n*m)))
end
figure, subplot(1,2,1), semilogy(s,'k')     % valores singulares
subplot(1,2,2), plot(cumsum(s)/sum(s),'k')  % soma acumulada

Guardar \(\mathbf{\tilde{U}}\), \(\boldsymbol{\tilde{\Sigma}}\) e \(\mathbf{\tilde{V}}\) exige \(r(n + m + 1)\) números, contra \(nm\) da imagem original: a truncação é uma compressão.

Imagem reconstruída para vários postos \(r\)

Com \(r = 5\) só se reconhecem as grandes regiões claras e escuras; com \(r = 100\) a imagem já é visualmente próxima da original.

Valores singulares e erro de truncamento

r =   5:  soma acumulada = 42.4%,  erro relativo (1.8) = 2.62%
r =  20:  soma acumulada = 57.0%,  erro relativo (1.8) = 0.77%
r = 100:  soma acumulada = 75.7%,  erro relativo (1.8) = 0.14%