Rank-one Riemannian Subspace Descent for Nonlinear Matrix Equations

Enunciado do Problema
O artigo aborda o desafio computacional de encontrar soluções de Matrizes Positivas Definidas Simétricas (SPD) para equações matriciais não lineares (NMEs) que surgem na teoria de controle, programação dinâmica e filtragem estocástica. Estas equações, que incluem as Equações Algébricas de Riccati de Tempo Contínuo e Discreto (CARE/DARE) e equações de Lyapunov, são frequentemente formuladas como a busca por $X \in \mathcal{P}_n$ tal que $G(X, A, C, \dots) = 0$. Isto é reestruturado como um problema de mínimos quadrados não lineares: $\min_{X \in \mathcal{P}_n} \|G(X, A, C, \dots)\|_F^2$.

O principal gargalo é a escalabilidade. Solvers padrão para equações matriciais densas tipicamente incorrem em $O(n^3)$ operações de ponto flutuante (flops) por iteração. Embora existam métodos $O(n^2)$, eles geralmente dependem de suposições de esparsidade ou soluções de baixo posto, que não se aplicam a problemas densos gerais. Consequentemente, para dimensões de estado elevadas (por exemplo, $n=10.000$), os solvers densos padrão tornam-se computacionalmente inviáveis devido a restrições de tempo e memória.

Metodologia: Descida de Subespaço Riemanniano de Posto Um (R1RSD)
Os autores propõem um algoritmo de Descida de Subespaço Riemanniano de Posto Um (R1RSD) que opera diretamente na variedade das matrizes SPD ($\mathcal{P}_n$) equipada com uma métrica riemanniana afim-invariante.

1. Seleção de Subespaço via Decomposição Espectral: Diferente dos métodos tradicionais de descida de subespaço que projetam gradientes em subespaços de coordenadas fixos, o R1RSD decompõe o gradiente riemanniano transformado em uma soma de componentes autótonas de posto um. O algoritmo seleciona o componente autótono dominante (correspondente ao maior valor absoluto do autovalor) para definir a direção de descida.
2. Método de Potência para Identificação de Direção: Para evitar o custo $O(n^3)$ de uma decomposição de autovalores completa, o componente dominante é identificado usando o método da potência. Os autores demonstram que um pequeno número de iterações de potência (especificamente $r = O(\log n)$) é suficiente para aproximar efetivamente a direção dominante.
3. Atualizações de Posto Um e Complexidade: A regra de atualização é derivada usando a exponencial de matriz de uma matriz de posto um, que admite uma atualização de posto um em forma fechada. Ao manter o fator de Cholesky $B_t$ de $X_t$ e seu inverso via atualizações de posto um (usando a identidade de Woodbury e cholupdate), o algoritmo evita multiplicações de matrizes densas.
4. Produtos Gradiente-Vetor: Para objetivos definidos como composições de operações matriciais padrão (adição, multiplicação, inversa, log-determinante), o produto gradiente-vetor $\nabla f(X)v$ é computado em $O(n^2)$ de tempo através do travessia do grafo computacional e aplicação de produtos matriz-vetor, evitando a formação explícita da matriz gradiente completa.
5. Seleção de Tamanho de Passo: A estrutura de posto um permite a busca de linha exata em muitos casos com custo adicional mínimo, embora tamanhos de passo fixos ou de Armijo sejam viáveis.

Principais Contribuições

• Inovação Algorítmica: O artigo introduz um método que atualiza as iterações ao longo de uma única direção de posto um cuidadosamente selecionada derivada do gradiente riemanniano, diferenciando-se de métodos que projetam em subespaços fixos.
• Redução de Complexidade: O algoritmo proposto alcança um custo aritmético de $O(n^2 \log n)$ por iteração (dominado por $r$ iterações de potência de $O(n^2)$ cada) e requer apenas $O(n)$ iterações sob suposições padrão. Isto contrasta com métodos de descida de subespaço de estado da arte (ex: [43]) que possuem custo de $O(n^2)$ por iteração mas requerem $O(n^2)$ iterações totais, e solvers padrão que requerem $O(n^3)$ por iteração.
• Garantias Teóricas: Os autores fornecem análise de convergência tanto para casos geo-suaves (não convexos) quanto para casos geo-fortemente convexos. Eles estabelecem um limite de $O(n)$ iterações para o caso geral e um limite melhorado para o caso fortemente convexo, com contagens totais de flops de $O(n^3 \log n)$ e $O(n^3 \kappa \log n)$, respectivamente.
• Implementação Prática: O trabalho inclui uma implementação MATLAB disponível publicamente que utiliza atualizações de baixo posto para manter a eficiência, garantindo que a complexidade de $O(n^2)$ por iteração seja realizada na prática.

Resultados Numéricos
Os autores comparam o R1RSD com solvers de estado da arte, incluindo o icare do MATLAB, algoritmos de Dobramento de Preservação de Estrutura (SDA) e o método de descida de coordenada baseado em Fatoração Burer-Monteiro (BMFC).

• Escalabilidade: Experimentos em CARE, DARE e outras NMEs com dimensões de até $n=10.000$ demonstram que os solvers padrão (SDA, icare) e o BMFC falham em retornar soluções ou esgotam a memória para $n$ elevado. O R1RSD resolve com sucesso estas instâncias.
• Eficiência: Para $n=100$, o R1RSD converge em significativamente menos iterações (1 a 2 ordens de magnitude) comparado ao BMFC.
• Robustez: O algoritmo permanece eficaz mesmo quando o tamanho de passo ótimo varia significativamente, e a seleção gulosa de subespaço (usando $\approx 10$ iterações de potência) fornece um forte equilíbrio entre velocidade de convergência e custo por iteração.

Significância e Alegações
O artigo afirma que o R1RSD provê uma abordagem prática para equações matriciais densas de alto nível e com restrição SPD onde solvers densos tradicionais são proibitivos. Ao reduzir a complexidade por iteração para $O(n^2)$ enquanto mantém um número de iterações razoável, o método preenche a lacuna entre a escalabilidade teórica e a aplicabilidade prática para problemas de controle de grande escala. Os autores observam que esta abordagem é particularmente relevante para cenários de controle de tempo variável e adaptativo (como estruturas de Equação de Riccati Dependente do Estado) onde resolver equações de Riccati a cada passo de tempo seria, de outra forma, inviável. O trabalho enfatiza que as atualizações de variedade de posto um podem lidar efetivamente com problemas que os solvers padrão não conseguem, sem depender de suposições de esparsidade ou de baixo posto sobre as soluções ou coeficientes.