DMD — Algoritmo Exato

De Onde Vem o \(K\) da Introdução?

Em Koopman — Introdução, \(A\) e \(K\) apareciam como soluções de um mínimos quadrados, sem mostrar como calculá-las na prática quando \(n\) (nº de estados/sensores) é grande. O DMD exato (Tu et al. 2014) resolve exatamente esse problema — é o caso \(\phi(x)=x\) do EDMD, e a base algorítmica de tudo que vem depois.

• Formar \(A\in\mathbb{R}^{n\times n}\) diretamente é caro (ou nem cabe em memória) quando \(n\) é grande e o nº de amostras \(m\) é pequeno;
• Solução: nunca formar \(A\) — projetar num subespaço de posto baixo via SVD e recuperar os autovetores de \(A\) exatamente, sem aproximação.

Formalismo: Definindo o Operador DMD

Dado um par de sequências de instantâneos (snapshots) com um passo de avanço, \(x_k'=f(x_k)\), montam-se as matrizes

\[ X = \begin{bmatrix}x_1 & x_2 & \cdots & x_m\end{bmatrix} \in\mathbb{R}^{n\times m}, \qquad X' = \begin{bmatrix}x_1' & x_2' & \cdots & x_m'\end{bmatrix} \in\mathbb{R}^{n\times m} \]

O operador DMD é o \(A\) que melhor explica a transição linearmente, no sentido de mínimos quadrados:

\[ A = \arg\min_{A} \|X' - AX\|_F = X'X^{+} \]

onde \(X^{+}\) é a pseudo-inversa de Moore–Penrose de \(X\).

Formalismo: O Obstáculo da Dimensão

\(A\in\mathbb{R}^{n\times n}\) — para \(n\) na casa de centenas ou milhares (sensores, malha espacial, estados de um conversor multi-nível), formar, guardar e autodecompor \(A\) é inviável, mesmo com \(m\ll n\) amostras.

Formalismo: Projetando o Operador (\(\tilde A\))

Substituindo \(X^{+}\) na definição de \(A\): \(A = X'V\Sigma^{-1}U^{*}\). Define-se o operador projetado na base POD \(U\) como \(\tilde A = U^{*}AU\). Expandindo passo a passo:

\[ \tilde A = U^{*}\big(X'V\Sigma^{-1}U^{*}\big)U = U^{*}X'V\Sigma^{-1}\big(U^{*}U\big) = U^{*}X'V\Sigma^{-1} \]

usando \(U^{*}U=I_r\). Logo

\[ \boxed{\tilde A = U^{*}X'V\Sigma^{-1}} \in\mathbb{C}^{r\times r} \]

barato de autodecompor mesmo quando \(A\) (\(n\times n\)) nunca chega a ser formado.

Formalismo: Reconstrução Exata dos Modos

Prova, substituindo \(A=X'V\Sigma^{-1}U^{*}\):

\[ A\varphi_i = \big(X'V\Sigma^{-1}U^{*}\big)\big(X'V\Sigma^{-1}\big)w_i = X'V\Sigma^{-1}\underbrace{\big(U^{*}X'V\Sigma^{-1}\big)}_{\tilde A}w_i = X'V\Sigma^{-1}\lambda_i w_i = \lambda_i\varphi_i \]

Os \(\varphi_i\) são os modos DMD exatos — diferem dos modos projetados \(Uw_i\), que só coincidem com autovetores de \(A\) quando as colunas de \(X'\) pertencem ao espaço-coluna de \(U\).

Algoritmo — Pseudocódigo

Algoritmo 1 — DMD Exato (Tu, Rowley, Luchtenburg, Brunton & Kutz, 2014)
Entrada: instantâneos x_1,...,x_{m+1} ∈ R^n; passo Δt; posto r
Saída:   modos Φ, autovalores Λ, amplitudes b

 1: X  ← [x_1  x_2  ...  x_m]
 2: X' ← [x_2  x_3  ...  x_{m+1}]
 3: (U, Σ, V) ← SVD(X), truncada aos r maiores valores singulares
 4: Ã ← U* X' V Σ⁻¹                        # operador projetado, r×r
 5: (W, Λ) ← autodecomposição(Ã)           # Ã W = W Λ
 6: Φ ← X' V Σ⁻¹ W                          # modos DMD exatos
 7: para cada i: ω_i ← ln(λ_i) / Δt        # taxa de crescimento contínua
 8: b ← Φ⁺ x_1                              # amplitudes iniciais
 9: retornar Φ, Λ, b

Algoritmo — Escolha do Posto \(r\)

• Os valores singulares \(\sigma_i\) (linha 3) decaem — escolher \(r\) pelo “cotovelo” ou pela energia acumulada \(\sum_{i\le r}\sigma_i^2\big/\sum_i\sigma_i^2 \ge\) limiar (ex.: 99,9%);
• \(r\) pequeno demais descarta dinâmica relevante; \(r\) grande demais ajusta ruído — \(\tilde A\) fica mal-condicionada;
• No exemplo numérico a seguir os dados têm posto exato 2: o corte é imediato, sem ambiguidade.

Algoritmo — Modos Exatos vs. Projetados

Exemplo: Oscilador Harmônico Amortecido

\[ \ddot x + 2\zeta\omega_n\dot x + \omega_n^2 x = 0, \qquad x(t) = e^{-\sigma t}\cos(\omega_d t), \quad \sigma=\zeta\omega_n,\;\; \omega_d=\omega_n\sqrt{1-\zeta^2} \]

Estado \(s(t)=[x(t),\,\dot x(t)]\) é observado por \(n=30\) “sensores” virtuais — combinações lineares aleatórias de \(s(t)\), como amostras espacialmente distribuídas de um único modo físico:

import numpy as np

zeta, wn, dt = 0.05, 2*np.pi*1.0, 0.02
sigma, wd = zeta*wn, wn*np.sqrt(1 - zeta**2)

m = 120
t = np.arange(m + 1) * dt
x = np.exp(-sigma*t) * np.cos(wd*t)
v = np.exp(-sigma*t) * (-sigma*np.cos(wd*t) - wd*np.sin(wd*t))
S = np.vstack([x, v])                       # 2 x (m+1)

n = 30
rng = np.random.default_rng(0)
C = rng.normal(size=(n, 2))
Y = C @ S                                    # n x (m+1), posto exato 2

X, Xp = Y[:, :-1], Y[:, 1:]                  # linhas 1-2 do Algoritmo 1

Exemplo: Aplicando o Algoritmo (linhas 3–8)

U, Sig, Vt = np.linalg.svd(X, full_matrices=False)
r = 2
Ur, Sigr, Vr = U[:, :r], Sig[:r], Vt[:r, :].conj().T

Atil = Ur.conj().T @ Xp @ Vr @ np.diag(1/Sigr)      # linha 4
lam, W = np.linalg.eig(Atil)                        # linha 5
Phi = Xp @ Vr @ np.diag(1/Sigr) @ W                 # linha 6

omega = np.log(lam) / dt                            # linha 7
b = np.linalg.lstsq(Phi, Y[:, 0], rcond=None)[0]     # linha 8

print("omega estimado:", np.round(omega, 4))

omega estimado: [-0.3142+6.2753j -0.3142-6.2753j]

print("(-sigma, omega_d) verdadeiro:", (round(-sigma, 4), round(wd, 4)))

(-sigma, omega_d) verdadeiro: (-0.3142, np.float64(6.2753))

Exemplo: Modos e Reconstrução

import matplotlib.pyplot as plt

recon = np.real(Phi @ (b[:, None] * lam[:, None]**np.arange(m + 1)))
erro = np.max(np.abs(recon - Y))
print("erro máx. de reconstrução:", erro)

erro máx. de reconstrução: 1.5898393712632242e-13

fig, (ax1, ax2) = plt.subplots(1, 2)

theta = np.linspace(0, 2*np.pi, 200)
_ = ax1.plot(np.cos(theta), np.sin(theta), "k--", lw=1)
_ = ax1.plot(lam.real, lam.imag, "o", ms=10, color="#96263a")
_ = ax1.set_xlabel(r"Re($\lambda$)"); ax1.set_ylabel(r"Im($\lambda$)")
_ = ax1.set_title("Autovalores DMD"); ax1.set_aspect("equal")

_ = ax2.stem(np.arange(n), Phi[:, 0].real)
_ = ax2.set_xlabel("índice do sensor"); ax2.set_ylabel(r"Re($\varphi_1$)")
_ = ax2.set_title("Modo DMD dominante")

plt.tight_layout(); plt.show()

Ponte: Do DMD à Teoria de Koopman

DMD exato é o EDMD com o dicionário mais simples possível, \(\phi(x)=x\) — por isso só captura dinâmica já linear (ou linearizável) nas coordenadas medidas, como o oscilador acima.