← Matemática

Ortogonalização de Gram-Schmidt

Escolha a proximidade entre as colunas e observe as duas grelhas. Elas indicam o ângulo entre cada par de vetores devolvido pelas duas versões.

A carregar a simulação interativa...

Gram-Schmidt clássico devolve dois vetores separados por 60°

Prima ε = 10⁻⁸, clássico. Na grelha clássica, duas células indicam 60,00° para o par q₂ e q₃ nos dois cantos da matriz. Qualquer outro ângulo em ambas as grelhas marca 90,00°. O produto interno subjacente é 0,5000: um meio, com o rigor de uma unidade na última casa de um double.

A matriz tem três colunas, sendo cada uma delas o vetor (1, 0, 0, 0) com um único 10⁻⁸ inserido numa posição diferente. Peter Läuchli utilizou-a em 1961 e, desde então, ostenta o seu nome. Colunas quase repetidas constituem o caso problemático para qualquer ajuste de mínimos quadrados. Elas surgem sempre que duas das suas medições são praticamente idênticas.

O Gram-Schmidt modificado recebe as mesmas três colunas e todos os ângulos que devolve medem 90,00°, com o mesmo número de multiplicações.

Recue a separação um nível, para 10⁻⁷, e a mesma célula passará a marcar 89,24°. Um clique no controlo deslizante é a diferença entre errar por três quartos de grau ou errar por trinta.

As duas versões interpretam um número de forma diferente

Ajuste o guia para a coluna 3 e prima Clássico, seguido de Modificado. A linha O que q₂ subtrai à coluna 3 passa de 0 para 7,0711 × 10⁻⁹. Com três colunas, este é o único coeficiente que os dois algoritmos calculam de forma distinta. A entrada diagonal por baixo dele, tal como todos os ângulos da grelha, resultam dessa diferença.

A versão clássica avalia esse coeficiente como q₂ᵀa₃, em relação à coluna original que recebeu. A versão modificada avalia q₂ᵀv, em relação ao que sobra da coluna após a remoção da parte correspondente a q₁. No papel, tratam-se do mesmo número. O vetor v difere de a₃ apenas por um múltiplo de q₁, e q₂ é perpendicular a q₁. Na aritmética de máquina, q₂ não é perfeitamente perpendicular a q₁, e a₃ ainda preserva 7,0711 × 10⁻⁹ de si próprio na direção de q₂. A versão clássica nunca inspeciona esse aspeto, logo q₃ retém esse valor, o que origina os 60°.

Arraste o controlo de 10⁻¹ até 10⁻⁸ e observe a separação das duas curvas no painel inferior. Ajustada sobre 113 pontos nesse intervalo, a perda de ortogonalidade da versão clássica cresce como κ1,98, enquanto a da versão modificada cresce como κ0,98. Uma simples reordenação na subtração vale toda uma potência do número de condição.

A base fica arruinada, mas a fatorização é exata

A última linha da tabela de resultados é ‖A − QR‖ / ‖A‖. Para ε = 10⁻⁸, com a versão Clássica, lê-se 0. Multiplique o Q arruinado pelo R que o acompanha, e a matriz original é reconstituída, entrada a entrada, sem qualquer resíduo.

Percorrendo todos os 280 estados possíveis nos controlos — quatro contagens de colunas, dez separações, dois algoritmos e todas as colunas de cada um —, esse resíduo relativo nunca ultrapassa 5,5511 × 10⁻¹⁷. O intervalo entre 1 e o double seguinte é de 2,2204 × 10⁻¹⁶. O valor mantém-se, portanto, num quarto de uma unidade na última casa.

O método de Gram-Schmidt efetivamente escreveu A como um produto dos vetores que devolveu pelos coeficientes que registou. O que deixou de fazer foi devolver vetores perpendiculares entre si. Saber qual destas duas propriedades é a necessária será uma questão pertinente para o seu próprio código. A projeção de um vetor com Qᵀb apoia-se na propriedade que falhou, enquanto a reconstrução de A a partir de Q e R apoia-se na que se manteve.

Problema resolvido na íntegra

  1. Onde o Gram-Schmidt clássico perde o ângulo reto 8 passos

    A ferramenta apresenta três colunas separadas por 10⁻⁸ e a respetiva resposta do Gram-Schmidt clássico. Execute o algoritmo à mão, encontre o arredondamento isolado que lhe custa trinta graus e, em seguida, indique até onde cada versão pode ser forçada antes de deixar de valer a pena executá-la.

    Abrir este problema na ferramenta

    1. Cada coluna corresponde ao vetor (1, 0, 0, 0) com um ε colocado numa posição diferente. A matriz AᵀA é a matriz de uns de dimensões 3 × 3 somada a ε²I. Os seus valores próprios são 3 + ε² uma vez e ε² duas vezes. Como tal, os valores singulares são √(3 + ε²) e ε, o que faz com que o número de condição seja exato e não uma estimativa. A primeira linha da tabela regista 1,7321 × 10⁸.

      A=(111ε000ε000ε),ε=10−8,κ(A)=3+ε2ε=1.7321×108A=\begin{pmatrix}1&1&1\\\\ \varepsilon&0&0\\\\ 0&\varepsilon&0\\\\ 0&0&\varepsilon\end{pmatrix}, \qquad \varepsilon=10^{-8}, \qquad \kappa(A)=\frac{\sqrt{3+\varepsilon^{2}}}{\varepsilon}=1{.}7321\times10^{8}

    2. ‖a₁‖² é 1 + ε², sendo ε² igual a 10⁻¹⁶. Um double armazena o 1, e o próximo número representável encontra-se 2,2204 × 10⁻¹⁶ mais à frente. Qualquer valor inferior a metade dessa distância, 1,1102 × 10⁻¹⁶, desaparece ao ser somado a 1. Dez elevado a menos dezasseis fica abaixo desse limiar. A raiz quadrada é 1, dividir por este valor não altera nada, pelo que q₁ surge como a própria coluna, com o ε² eliminado do registo. Todos os outros arredondamentos neste problema custam uma unidade na última casa. Este custa trinta graus.

      ∥a1∥2=1+ε2=1+10−16,10−16<1.1102×10−16 ⟹ fl(1+ε2)=1,q1=a1\lVert a_{1}\rVert^{2}=1+\varepsilon^{2}=1+10^{-16}, \qquad 10^{-16}<1{.}1102\times10^{-16} \ \Longrightarrow\ \mathrm{fl}\bigl(1+\varepsilon^{2}\bigr)=1, \qquad q_{1}=a_{1}

    3. Passando à coluna 2, ambas as versões concordam, dado existir apenas um vetor para trás. O coeficiente de projeção é q₁ᵀa₂ = 1. A sua subtração deixa (0, −ε, ε, 0), que possui um comprimento de ε√2 = 1,4142 × 10⁻⁸. A coluna 2 chegou com comprimento 1 até onde um double consegue distinguir, e quase toda ela já estava coberta por q₁, à exceção de 1,4142 × 10⁻⁸.

      r12=q1Ta2=1,v2=a2−q1=(0, −ε, ε, 0)T,∥v2∥=ε2=1.4142×10−8r_{12}=q_{1}^{\mathsf{T}}a_{2}=1,\qquad v_{2}=a_{2}-q_{1}=(0,\,-\varepsilon,\,\varepsilon,\,0)^{\mathsf{T}},\qquad \lVert v_{2}\rVert=\varepsilon\sqrt{2}=1{.}4142\times10^{-8}

    4. Coluna 3, classicamente. O coeficiente em relação a q₁ volta a ser 1. O coeficiente em relação a q₂ é q₂ᵀa₃. Visto que q₂ existe apenas nas posições 2 e 3, enquanto a₃ mantém o seu ε na posição 4, este produto interno é 0. É precisamente isso que indica a linha O que q₂ subtrai à coluna 3. Subtraindo apenas q₁, resta (0, −ε, 0, ε), novamente com o comprimento ε√2. A linha O que resta da coluna 3 também regista 1,4142 × 10⁻⁸.

      r13=q1Ta3=1,r23=q2Ta3=0,v3=a3−q1=(0, −ε, 0, ε)T,q3=12(0, −1, 0, 1)Tr_{13}=q_{1}^{\mathsf{T}}a_{3}=1,\qquad r_{23}=q_{2}^{\mathsf{T}}a_{3}=0,\qquad v_{3}=a_{3}-q_{1}=(0,\,-\varepsilon,\,0,\,\varepsilon)^{\mathsf{T}},\qquad q_{3}=\tfrac{1}{\sqrt{2}}(0,\,-1,\,0,\,1)^{\mathsf{T}}

    5. Tanto q₂ como q₃ foram construídos subtraindo q₁ a uma coluna que partilha o seu 1 inicial. Consequentemente, ambos transportam −1/√2 na posição 2 e não partilham mais nada. O produto interno entre eles é um meio. A ferramenta imprime 0,5000 e 60,00°, mas o valor que armazena é 0,4999999999999999, que dista uma unidade na última casa do meio exato garantido pela álgebra. Isto acontece porque ±1/√2 não constitui uma fração binária. Esta é a magnitude de erro que tal aritmética consegue cometer. Os trinta graus não advêm daqui. A linha de perda da versão clássica regista 0,7071, o que equivale a √(2 × 0,5²): o mesmo meio, contabilizado uma vez em cada canto de QᵀQ.

      q2Tq3=12[(−1)(−1)+(1)(0)]=12,arccos⁡12=60∘,2⋅0.52=0.7071q_{2}^{\mathsf{T}}q_{3} =\tfrac{1}{2}\bigl[(-1)(-1)+(1)(0)\bigr]=\tfrac{1}{2},\qquad \arccos\tfrac{1}{2}=60^{\circ},\qquad \sqrt{2\cdot 0{.}5^{2}}=0{.}7071

    6. Na versão modificada, para a mesma coluna. Após remover a parte de q₁, o resto é (0, −ε, 0, ε), e esta versão faz as projeções contra este vetor em vez de a₃. O produto interno de q₂ com ele resulta em ε/√2 = 7,0711 × 10⁻⁹, e não zero. Subtraindo este valor, o que sobra é (0, −ε/2, −ε/2, ε), com comprimento ε√(3/2) = 1,2247 × 10⁻⁸. O cálculo de q₂ᵀq₃ dá agora 1,1102 × 10⁻¹⁶ em vez de um meio. A linha de perda modificada regista 1,1547 × 10⁻⁸. Tudo o que aqui resta é o que q₁ não conseguiu retirar: a subtração deixa −ε² na direção de q₁, e a divisão por um resto com comprimento de apenas ε√2 volta a transformar isso num ε/√2.

      r23=q2Tv3=ε2=7.0711×10−9,v3↦(0, −ε2, −ε2, ε)T,∥v3∥=ε32=1.2247×10−8,q2Tq3=0r_{23}=q_{2}^{\mathsf{T}}v_{3}=\frac{\varepsilon}{\sqrt{2}}=7{.}0711\times10^{-9},\qquad v_{3}\mapsto\bigl(0,\,-\tfrac{\varepsilon}{2},\,-\tfrac{\varepsilon}{2},\,\varepsilon\bigr)^{\mathsf{T}},\qquad \lVert v_{3}\rVert=\varepsilon\sqrt{\tfrac{3}{2}}=1{.}2247\times10^{-8},\qquad q_{2}^{\mathsf{T}}q_{3}=0

    7. O mesmo raciocínio aplica-se a qualquer contagem de colunas disponibilizada pela ferramenta. A versão clássica constrói cada q subsequente subtraindo q₁ a uma coluna com a qual partilha o seu 1 inicial. Portanto, todos os pares entre eles saem a 60°: são (n − 1)(n − 2)/2 pares, com uma perda de ortogonalidade igual à raiz quadrada de metade dessa contagem. Com duas colunas não existem tais pares. O que sobra é a inclinação da ordem de ε partilhada por ambas as versões.

      nκ(A)#{60∘}∥QCTQC−I∥F∥QMTQM−I∥F21.4142×10801.0000×10−81.0000×10−831.7321×10810.70711.1547×10−842.0000×10831.22471.2247×10−852.2361×10861.73211.2649×10−8\begin{array}{r|rrrr} n & \kappa(A) & \#\{60^{\circ}\} & \lVert Q_{C}^{\mathsf{T}}Q_{C}-I\rVert_{F} & \lVert Q_{M}^{\mathsf{T}}Q_{M}-I\rVert_{F} \\\\ \hline 2 & 1{.}4142\times10^{8} & 0 & 1{.}0000\times10^{-8} & 1{.}0000\times10^{-8} \\\\ 3 & 1{.}7321\times10^{8} & 1 & 0{.}7071 & 1{.}1547\times10^{-8} \\\\ 4 & 2{.}0000\times10^{8} & 3 & 1{.}2247 & 1{.}2247\times10^{-8} \\\\ 5 & 2{.}2361\times10^{8} & 6 & 1{.}7321 & 1{.}2649\times10^{-8} \end{array}

    8. A reta ajustada para a versão clássica atinge uma perda de 1 em κ = 2,35 × 10⁸, motivo pelo qual a curva clássica já se encontra plana na extremidade direita do painel. A versão modificada obtém um limite medido em vez de extrapolado, uma vez que a sua curva vira para lá de 10⁻⁸. Uma reta ajustada prolongada até esse ponto seria refutada pelo painel que se lhe segue. Ao longo de 724 pontos do varrimento, a sua perda nunca excedeu 1,25 épsilons de máquina vezes κ, e este limite só alcança o valor 1 em κ = 3,60 × 10¹⁵. O inverso do épsilon de máquina é 4,50 × 10¹⁵.

      log⁡10κ=16.57551.9801=8.3710 ⟹ κ=2.35×1081.25⋅(2.2204×10−16) κ=1 ⟹ κ=3.60×1015(2.2204×10−16)−1=4.50×1015\begin{array}{rcl} \log_{10}\kappa=\dfrac{16{.}5755}{1{.}9801} &=& 8{.}3710 \ \Longrightarrow\ \kappa=2{.}35\times10^{8} \\\\ 1{.}25\cdot\bigl(2{.}2204\times10^{-16}\bigr)\,\kappa &=& 1 \ \Longrightarrow\ \kappa=3{.}60\times10^{15} \\\\ \bigl(2{.}2204\times10^{-16}\bigr)^{-1} &=& 4{.}50\times10^{15} \end{array}

    Resposta

    Os trinta graus emergem de uma única linha. A versão clássica lê q₂ᵀa₃ = 0, enquanto a versão modificada lê q₂ᵀv₃ = 7,0711 × 10⁻⁹. No papel, trata-se do mesmo número.

    A separação de 10⁻⁸ corresponde ao primeiro nível em que ε² desaparece de 1 + ε², dado que √(metade do épsilon de máquina) é 1,0537 × 10⁻⁸. A partir daqui para baixo, a célula exibe 60,00° independentemente do quanto continue a arrastar. Um nível atrás, a 10⁻⁷, a leitura é de 89,24° e a versão clássica já perdeu 1,88 × 10⁻² de ortogonalidade. O colapso começa muito antes de o controlo deslizante chegar à marca que o torna exato.

    Åke Björck publicou a análise da versão modificada em 1967 e, desde então, a reordenação passou a ser a norma. Ambas executam o mesmo número de multiplicações, mas para esta separação a perda de ortogonalidade difere por um fator de 6,1 × 10⁷. A versão modificada mantém os seus ângulos retos até quase ao limite em que um double deixa de conseguir distinguir as colunas. A versão clássica cede sete ordens de grandeza mais cedo.

    q2Ta3=0≠q2Tv3=ε2=7.0711×10−9q_{2}^{\mathsf{T}}a_{3}=0 \qquad\neq\qquad q_{2}^{\mathsf{T}}v_{3}=\frac{\varepsilon}{\sqrt{2}}=7{.}0711\times10^{-9}
Referências (5)

· Verificado com fontes publicadas. Correções bem-vindas.