Decomposição QR

$\newcommand{\R}{\mathbb{R}}$

Na aula passada vimos que o problema de ajuste de curvas recai em encontrar o vetor $x$ para o qual a norma euclidiana do resíduo de um sistema linear é a menor possível. Esta é a maneira de definir a solução de quadrados mínimos para um sistema linear.

Com essa constatação, um sistema linear $Ax=b,$ com $A\in\R^{m\times n},$ $m\gt n$ (sistema sobredeterminado) tem como solução de quadrados mínimos o vetor $x^*$ que é solução do sistema normal $$A^TAx=A^Tb.$$

Nesta aula vamos ver outra forma de determinar esse vetor $x^*,$ que dispensa a construção e resolução do sistema normal. Para tanto, faremos uso da decomposição QR, cuja definição e existência é assegurada pelo teorema a seguir.

Teorema: Toda matriz $A \in \R^{m\times n}$ pode ser decomposta no produto de uma matriz ortogonal $Q\in\R^{m\times m}$ por uma matriz triangular superior $R\in\R^{m\times n},$ ou seja, $$ A = QR. $$

Mas o que é uma matriz ortogonal? Um matriz $Q,$ quadrada, é dita ortogonal se satisfaz a seguinte equação $$Q^TQ = I.$$ Desta definição, imediatamente surge a propriedade de que $Q^T$ é de fato a inversa de $Q.$ Outra propriedade particularmente importante para esta aula é que transformações por matrizes ortogonais não alteram a norma euclidiana (ou norma-2) de vetores. Com efeito, veja que $$ \|Qx\|_2^2 = (Qx)^T(Qx) = x^TQ^TQx=x^Tx= \|x\|_2^2, $$ onde usamos que $\|u\|_2^2 = u_1^2+\cdots+u_m^2 = u^Tu.$ Como $Q^T = Q^{-1},$ $QQ^T$ também é a matriz identidade. Com isso, é simples ver também que $\|Q^Tx\|_2^2 = \|x\|_2^2.$ Em resumo, temos as propriedades:

1. $Q^TQ=I$
2. $\|Qu\|_2^2 = \|u\|_2^2 = \|Q^Tu\|_2^2$

Mas como isso nos ajuda no problema de quadrados mínimos? Se $A=QR,$ observe que o resíduo que queremos minimizar pode ser escrito como \begin{align*} \|Ax-b\|_2^2 & = \|QRx-b\|_2^2 \\ & =\|Q^T(QRx-b)\|_2^2 &\mbox{(propriedade 2)}\\ & =\|Rx-Q^Tb\|_2^2 &\mbox{(propriedade 1)} \end{align*}

Como a matriz $R$ é triangular superior e tem as mesmas dimensões da matriz $A,$ $R$ pode ser escrita como $$ R = \begin{bmatrix}\hat{R}\\\ \mathbf{0}\end{bmatrix}, $$ onde $\hat{R}\in\R^{n\times n}$ é uma matriz triangular superior e $\mathbf{0}$ é a matriz nula de dimensão $(m-n)\times n.$

Repartindo o vetor $Q^Tb =z$ da mesma forma, temos que $$ \|Rx-Q^Tb\|_2^2 = \|{\left[\begin{array}{c}\hat{R}x-\hat{z}\\\ -\tilde{z}\end{array}\right]}\|_2^2 = \|\hat{R}x-\hat{z}\|_2^2+\|\tilde{z}\|_2^2. $$

Logo $$\min\|Ax-b\|^2_2 \Longleftrightarrow \min\|Rx-z\|^2_2 \Longleftrightarrow \min\left\{\|\hat{R}x -\hat{z}\|^2_2 + \|\tilde{z}\|^2_2\right\}. $$ Acontece que na última expressão temos o resíduo do sistema linear $\hat{R}x=\hat{z},$ que é um sistema linear quadrado, com matriz de coeficientes triangular superior. Se $\hat{R}$ for não singular, estes sistema tem solução, que é facilmente computada (veja esta aula). Ao resolver este sistema exatamente, conseguimos o menor resíduo possível — zero — e portanto temos a solução de norma euclidiana mínima para o problema original. De brinde, podemos ver também que o resíduo do sistema $Ax=b$ na solução de quadrados mínimos é $\|\tilde{z}\|_2.$

Para que isto funcione basta que $\hat{R}$ seja não singular, o que é equivalente a pedir que as colunas da matriz $A$ sejam linearmente independentes. Se isto não ocorrer, ainda conseguiremos resolver o problema, mas a solução de quadrados mínimos não será única.

O caminho para encontrar a solução de quadrados mínimos de um sistema linear sobredeterminado está sumarizado nos passos a seguir.

Nesta aula não vimos como computar a decomposição QR, mas tão somente como utilizá-la. Foge ao escopo desse curso a discussão de métodos para a obtenção da decomposição QR. Para quem se interessar, recomendo os livros Greenbaum e Chartier (2012) e Steweart (1998).

No Octave, o comando [Q,R] = qr(A), retorna as matrizes $Q$ e $R$ da decomposição. De fato, o comando \, quando aplicado a um sistema sobredeterminado, realiza os três passos acima e retorna a solução no sentido de quadrados mínimos. Sendo assim, se tiver acesso ao Octave, pode usar diretamente o \ para resolver sistemas sobredeterminados, sabendo agora o que ele faz por debaixo dos bits.

Para finalizar, vejamos um exemplo no Octave.

A = [1 1; 1 2; 2 3];
b = [ 2; 3; 4];
[Q,R] = qr(A)
Q =

  -4.0825e-01   7.0711e-01  -5.7735e-01
  -4.0825e-01  -7.0711e-01  -5.7735e-01
  -8.1650e-01   5.5511e-16   5.7735e-01

R =

  -2.44949  -3.67423
   0.00000  -0.70711
   0.00000   0.00000

z = Q'*b
z =

  -5.30723
  -0.70711
  -0.57735

RR = R(1:2,1:2);
zz = z(1:2,1);
x = RR\zz
x =

   0.66667
   1.00000

A\b
ans =

   0.66667
   1.00000

norm(A*x-b)
ans =  0.57735

Como exercício, repita o bloco acima no Octave, exibindo as saídas intermediárias (retirando os ; do final de cada linha), para verificar se ficou claro o propósito de cada linha executada.

Referência

Anne Greenbaum e Timothy P. Chartier. Numerical Methods: Design, Analysis, and Computer Implementation of Algorithm. Princeton University Press, 2012.

Gilbert W. Steweart. Afternotes goes to Graduate School: Lectures on Advanced Numerical Analysis. SIAM, 1998.

1. Verique que matrizes abaixo são ortogonais.

1. $A = \displaystyle{\left[\begin{array}{rr}0 & -1\\\ -1 &0\end{array}\right]}$
2. $B = \displaystyle{{1\over{\sqrt{2}}}\left[\begin{array}{rr}-1 & 1\\\ 1 & 1\end{array}\right]}$
3. $R = \displaystyle{\left[\begin{array}{rr}\cos x & \sin x\\\ -\sin x & \cos x\end{array}\right]}$
4. $H = I-2uu^T$, onde $I$ é a matriz identidade e $u$ é um vetor tal que $u^Tu=1$.

2. Seja $C=A^TB$, onde $A$ e $B$ são matrizes com dimensões adequadas para que o produto esteja bem definido.

1. Verifique que $c_{ij} = a_i^Tb_j,$ onde $c_{ij}$ é o elemento na linha $i$ e coluna $j$ da matriz $C,$ $a_i$ representa a coluna $i$ da matriz $A$ e $b_j$ representa a coluna $j$ da matriz $B.$
2. Se $Q$ é uma matriz ortogonal, usando o item anterior, verifique que suas colunas são ortogonais entre si e que cada coluna tem norma euclidiana unitária.

Para esse exercício, lembre que o produto interno usual entre dois vetores do $\R^n$ é definido como $\langle x,y\rangle = \sum_k x_ky_k = x^Ty$ e que a norma euclidiana de um vetor é dada por $\|u\|_2 = \sqrt{u^Tu}.$

3. Considere os sistemas lineares a seguir. Para cada um deles, encontre a solução de quadrados mínimos usando a decomposição QR e depois verifique seu resultado comparando-o com o encontrado diretamente pelo \. Compute o resíduo em cada caso, sem calcular o produto $Ax^*.$

1. $\displaystyle{\begin{bmatrix}1&2\\ 2&1\\ 1&1\end{bmatrix} \begin{bmatrix}x_1\\x_2\end{bmatrix} = \begin{bmatrix}3 \\ 3\\ 1\end{bmatrix}}$
2. $\displaystyle{\left[\begin{array}{rr}1&3\\ -2&1\\ 1&1\end{array}\right] \begin{bmatrix}x_1\\x_2\end{bmatrix} = \left[\begin{array}{r}1 \\ 5\\ -1\end{array}\right]}$
3. $\displaystyle{\begin{bmatrix}1\\ 2\\ 3\end{bmatrix} x = \begin{bmatrix}1 \\ 1\\ 1\end{bmatrix}}$
4. $\displaystyle{\begin{bmatrix}1\\ 1\end{bmatrix} x = \begin{bmatrix}1 \\ 2\end{bmatrix}}$