Introdução à Álgebra Linear Computacional:

Fatorações Básicas
Prof. Alexandre Salles da Cunha e Profa. Ana Paula Couto

2024/2

1 Introdução
Nesta seção apresentamos algumas fatorações matriciais básicas. Independen-
temente da forma de A ∈ Rm×n (quadrada, quando m = n, ou retangular
esbelta se n < m ou larga se n > m), fatorar uma matriz consiste no processo
algorı́tmico que permite escrever A como produto de outras matrizes com al-
guma propriede ou mesmo topologia particular, mais convenientes para algum
propósito especı́fico.
Por topologia queremos dizer padrão de esparsidade, ou seja, a localização de
uma região da matriz onde são autorizados a estarem localizados seus elementos
não nulos. Para uma dada topologia, fora dessa região especı́fica, todos os
elementos da matriz devem ser nulos. Dois exemplos de topologias de matrizes
são matrizes triangulares inferiores e superiores, definidas como as matrizes
que tem zeros em todas suas entradas acima e abaixo da diagonal principal,
respectivamente. Algumas propriedades de interesse dos fatores podem ser: ter
posto completo, terem colunas ou linhas ortonormais, por exemplo.
Ao longo de todo o curso de ALC, vamos discutir diversas fatorações e algo-
ritmos para computá-las. São elas:

1. Fatoração P A = LU , onde L é triangular inferior, com a diagonal unitária,

U é uma triangular superior e P é uma matriz de permutação.

2. Fatoração de Cholesky A = RT R, onde R é triangular inferior, com a

diagonal positiva. A precisa ser SPD.
3. Fatoração Espectral A é simétrica e A = QΛQT , onde Λ é uma matriz
diagonal com os autovalores de A e Q é ortogonal, como os autovetores de
A em suas colunas.

4. Fatoração de Schur A = QT QT onde Q é ortogonal, T é uma triangular

superior. A matriz A não precisa ser simétrica nem diagonalizável. A
diagonal de T armazena os autovalores de A.
5. Fatoração A = QR, onde Q é uma matriz com colunas ortonormais e R é
uma triangular superior.
6. Decomposição em Valores Singulares (Singular Value Decomposition -
SVD): A = U ΣV T , A ∈ Rm×n , U ∈ Rm×n , V ∈ Rn×n , Σ ∈ Rn×n ,

1
 V T V = In , U T U = Im e a matriz Σ é uma matriz de zeros, exceto pe-
las r primeiras entradas de sua diagonal, que guarda valores positivos
σ1 ≥ σ2 ≥ · · · σr , sendo r o posto de A.

Nesta seção vamos nos concentrar nas duas primeiras: fatoração P A = LU

e de Cholesky. Assim, os algoritmos que vamos desenvolver fatoram ou de-
compõem uma matriz A, cujas linhas foram trocadas de ordem por P , na forma
P A = LU onde L é uma triangular superior e U uma trinangular superior. A fa-
toração de Cholesky é um caso particular de P A = LU , onde os fatores L = U T
e P = I. A Fatoração de Cholesky só se aplica para matrizes simétricas po-
sitivas definidas, podendo ser adaptada para fatoração de matrizes de posto
incompleto, desde que sejam simétricas positivas semidefinidas. Isto é, podem
ser adaptadas para se fatorar A = RT R onde R ∈ Rn×r é triangular inferior
com Rii > 0, i = 1, . . . , r. Nessa seção vamos tratar primordialmente o caso
em que A é quadrada, embora a fatoração P A = LU pode ser adaptada para
produzir fatores para uma matriz retangular.
Nas fatorações que desejamos computar (seja P A = LU ou de Cholesky),
as matrizes L, U devem possuir posto completo, exatamente o posto de A (esta
sim, pode ter posto incompleto). Ou seja, as fatorações devem revelar o posto
da matriz A e apresentar bases para C(A), C(AT ).
Nossa opção é por denominar as duas fatorações estudadas nesta seção como
básicas, pelas seguintes razões:

1. Os elementos algorı́tmicos que empregam são bastante simples,

2. As bases fornecidas para C(A), C(AT ) não são ortonormais,
3. Com estas fatorações somos capazes de resolver boa parte dos sistemas
lineares com os quais nos deparamos em aplicações, desde que sejam bem
condicionados.

A fatoração P A = LU e de Cholesky são adequadas para se resolver sistemas

lineares que não sejam mal-condicionados e produzem resultados satisfatórios
para tais sistemas. Informalmente, sistemas lineares bem condicionados são
definidos por matrizes de coeficientes que, no processo de fatoração, não ten-
dem a gerar fatores com grande acúmulo de erros numéricos. Erros numéricos
grandes nos fatores se traduzem em erros numéricos grandes, por exemplo, nas
soluções dos sistemas lineares onde as matrizes aparecem. Para fatorar matrizes
mal-condicionadas, estudaremos fatorações especı́ficas, como a fatoração QR ou
SVD, na segunda metade do curso de ALC.

2 Razões para se fatorar matrizes

Antes de apresentarmos a primeira das fatorações discutidas aqui, é pertinente
mencionar razões para que fatoremos matrizes. Essencialmente, a fatoração
revela informação sobre a matriz A. No caso de sistemas dinâmicos lineares,
a fatoração revela informações sobre o sistema fı́sico que é representado pela
matriz.
Algumas das razões mais importantes para se fatorar uma matriz A ∈ Rm×n
são:

2
 1. Resolver um ou vários sistemas lineares, possivelmente definidos pela mesma
matriz de coeficientes A.
2. Analisar a existência e unicidade das soluções de sistemas lineares.
3. Calcular o determinante de uma matriz quadrada.
4. Conhecer espaços vetoriais associados à matriz: C(A), C(AT ), N (A), N (AT ).
Eventualmente, podemos desejar que as bases para estes espaços sejam
ortonormais. Para tanto, as fatorações empregadas devem levar estes as-
pectos em consideração.
5. Obter o espectro de A, ou seja seus autovetores e, eventualmente, seus
autovetores.
6. Conhecer os valores singulares de A, assim como seus vetores singulares,
de fundamental utilidade para o item abaixo.
7. Aproximar matrizes com muitas colunas ou muitas linhas por matrizes de
posto baixo. Com isso podemos resolver problemas aplicados da Ciência
da Computação em Otimização, em Inteligência Artificial, em Processa-
mento de Imagens e de Sinais, apenas para citar algumas aplicações.
8. As fatorações de matrizes nos permitem reformular problemas de Ma-
temática Aplicada, de uma forma mais conveniente, desde que o pro-
blema seja representado ou aproximado por um sistema linear.

Vamos brevemente discutir a primeira destas aplicações. Vamos supor que

precisemos resolver vários sistemas lineares Ax = b quadrados de ordem n, que
diferem entre si apenas pelo vetor de termos independenes. Isto é, a matriz de
coeficientes no sistema linear sempre é A, porém cada novo sistema linear possui
um novo vetor b. Vamos supor que posto de A seja completo e que A tenha sido
fatorada P A = LU . Vamos supor também que disponhamos de um algoritmo
capaz de resolver um sistema linear triangular, seja ele triangular superior ou
inferior. Para o desenvolvimento a seguir, recorde-se que P −1 = P T , quando
P é uma matriz de permutação. Veja que se dispomos da fatoração P A = LU
podemos substituir LU em P A e escrever:

Ax = b
P Ax = P b
LU x = P b
L(U x) = P b
Ly = P b (1)
Ux = y (2)

No desenvolvimento acima, usamos o fato de que U x é uma quantidade

desconhecida. Chamamos esta quantidade de y e então resolvemos o sistema
linear (1) com o algoritmo que supomos dispor. De posse desta quantidade y,
resolvemos o sistema (2), no qual y é agora conhecido e define o vetor de termos
independentes e x é a solução do sistema linear original.
Veja que, desde que disponhamos da fatoração P A = LU e que sejamos
capazes de resolver sistemas lineares triangulares, transformamos o problema

3
de encontrar Ax = b, onde A é uma matriz quadrada sem nenhuma topologia
particular, no problema de resolver dois sistemas lineares, todos os dois definidos
por matrizes triangulares: primeiro (1) e depois (2).
Vamos mostrar ao longo desta seção que o custo computacional de se fatorar
P A = LU é O(n3 ) e o custo de se resolver um sistema linear onde a matriz
de coeficientes é triangular (inferior ou superior, não importa) é O(n2 ). Então
suponha agora, que precisemos resolver k << n sistemas lineares distintos, de-
finidos pela mesma A. Fatoramos a matriz uma vez, pagando o custo de O(n3 ).
Para resolver os k sistemas lineares, resolvemos dois sistemas triangulares (1) e
(2), para cada um. Então somamos k(O(n2 ) + O(n2 )) ao custo computacional.
Ao fim, o custo computacional total é O(n3 + 2kn2 ) que é O(n3 ) para 2k << n.
Com a discussão acima, mostramos que podemos usar as fatorações do tipo
P A = LU e A = RT T (Cholesky) para resolvermos sistemas lineares de forma
eficiente.
Na próxima seção, vamos apresentar o algoritmo que resolve sistmas linea-
res triangulares inferiores e superiores, ingrediente importante para resolvermos
sistemas lineares mais gerais.

2.1 Resolução de sistemas lineares triangulares

Vamos começar esta seção formalizando as definições de matrizes e sistemas
lineares triangulares.
1. Uma matriz A, quadrada de ordem n, é triangular inferior se todos ele-
mentos acima da diagonal principal são nulos: aij = 0 para todo i, j : 1 ≤
i < j ≤ n.
2. Naturalmente, A é triangular superior se AT é triangular inferior.

3. Um sistema linear Ax = b é triangular, inferior ou superior, se a matriz

de coeficientes A é triangular, inferior ou superior, respectivamente.
Veja dois exemplos de matrizes triangulares.
Exemplo
 1 U é triangular
 superior
 e L é triangular
 inferior.
3 2 1 0 2 0 0 0
 L =  −1 2 0 0 
 0 1 2 3   
U =
 0 0 −2 1   3 1 −1 0 
0 0 0 4 4 1 −3 3
Como resolverı́amos o sistema linear Lx = b, onde b é um vetor qualquer,
por exemplo b = (2 3 2 9)T e L é a matriz do exemplo acima ? Vamos explorar
a esparsidade da matriz de coeficientes L e reescrever Lx = b de forma mais
conveniente, de forma que o algoritmo fique evidente. Usando o fato de que

4
lij = 0 : j > i para uma linha i do sistema linear Lx = b, temos:
n
X
lij xj = bi i = 1, . . . , n
j=1
i
X n
X
lij xj + lij xj = bi i = 1, . . . , n
j=1 j=i+1
i−1
X
lij xj + lii xi = bi i = 1, . . . , n
j=1
Pi−1
bi − j=1 lij xj
xi = k = 1, . . . , n
lii
Portanto, se calcularmos x1 , x2 , . . . , xi−1 nesta ordem, usando a expressão
(3) para k = i
Pk−1
bk − j=1 lkj xj
xk = , (3)
lkk
podemos calcular xi com as entradas já calculadas anteriormente x1 , x2 , . . . , xi−1 ,
necessárias na expressão (3). Este algoritmo, conhecido como Algoritmo de
Substituições Sucessivas, é apresentado abaixo na Figura 1.

function [y] = SubsSucessivas(L,b,n)

for i=1:n
soma = 0.0;
for k = 1:i-1
soma = soma + L(i,k)*y(k)
end
if (L(i,i) <> 0.0)
y(i) = (b(i)-soma)/L(i,i);
else
printf(’Matriz L e singular \n’)
break
end
end
endfunction

Figura 1: Algoritmo de Substituições Sucessivas.

Observe que se Lii = 0 para algum i, a matriz L e A são singulares. Pos-

teriormente vamos discutir como tratar o caso singular. Aqui, caso isso ocorra,
o algoritmo acusa a singularidade da matriz e interrompe sua execução. Veja o
exemplo de aplicação do algoritmo.

Exemplo 2 --> L,b

L =
2. 0. 0. 0.
-1. 2. 0. 0.

5
 3. 1. -1. 0.
4. 1. -3. 3.
b =
2.
3.
2.
9.
--> [y] = SubsSucessivas(L,b,size(L,1));
--> y’
ans =
1. 2. 3. 4.
Vamos agora mostrar que a complexidade do algoritmo de Substituições
Sucessivas é O(n2 ). O trecho de interesse do algoritmo, isto é, aquele que define
sua complexidade é indicado na Figura (2).
for i=1:n
soma = 0.0;
for k = 1:i-1
soma = soma + L(i,k)*y(k)
end
y(i) = (b(i)-soma)/L(i,i)
end

Figura 2: Trecho de interesse do algoritmo de Substituições Sucessivas.

Para avaliar a função de complexidade do algoritmo, vamos contar as operações

aritméticas de ponto flutuante realizadas pelo algoritmo. Não consideramos as
operações de incremento e comparação das variáveis inteiras, necessárias para as
estruturas for ou while, por exemplo. Cada operação aritmética de ponto flutu-
Pn
ante (+,-,×,÷) tem o mesmo custo unitário, 1. Lembramos que i=1 i = n(n+1) 2 .
Para um determinado valor de i fixo, a instrução
soma = soma + L(i,k)*y(k)
Pi−1
é executada g(i) = ( k=1 1) vezes, dentro da estrutura de controle
for k = 1:i-1
Pi−1
para este valor de i fixo. Veja que esta quantidade g(i) = k=1 1 é uma função
de i e assim sendo, para cada valor de i distinto, teremos uma contribuição
distinta. Cada vez que a instrução for executada, são realizadas uma soma e
uma multiplicação, ou seja, incorremos em um custo de 2 operações.
Agora, como podemos escrever o número de vezes que a estrutura de controle
for k = 1:i-1
é chamada ? Veja que esta estrutura de controle ela dentro de uma estrutura
de controle mais externa
for i=1:n
que controla os valores admissı́veis de i. Se desejamos o custo computacional
total relativo à instrução

6
soma = soma + L(i,k)*y(k)
precisamos avaliar g(1) + g(2) + · · · + g(n). Isso porque a estrutura mais interna
será executada tantas vezes quantos forem os valores assumidos de i. Então
podemos escrever a função de complexidade f (n) do algoritmo como:

n
X
f (n) = 2 g(i)
i=1
n X
X i−1
=2 1
i=1 k=1

Para resolver um somatório como o acima, começamos a explicitar o resul-

tado dos somatórios mais internos, pois estes assumem que valores fixos para as
variáveis foram definidos nos somatórios anteriores. Então temos:
n X
X i−1
f (n) = 2 1
i=1 k=1
Xn
=2 (i − 1)
i=1
 
n(n + 1)
=2 −n
2
= n(n + 1) − 2n
= n(n − 1)

Portanto, a instrução que estudadmos adiciona n(n − 1) operações aritméticas

de ponto flutuante (flops) ao custo computacional do algoritmo.
O custo total do algoritmo deve levar em conta também o custo adicionado
pela instrução
y(i) = (b(i)-soma)/L(i,i)
Pn
que é i=1 2 = 2n. Portanto, o custo total de substituições Sucessivas é n(n −
1) + 2n = n(n + 1), ou seja, seu custo pertence à classe de complexidade O(n2 ).

Vamos agora discutir a resolução de sistemas lineares triangulares supe-

riores. Uma vez que já discutimos os sistemas lineares inferiores e o algoritmo
de substituições que o resolve detalhadamente, nossa exposição do caso triangu-
lar superior é mais breve. Tudo é análogo ao caso triangular inferior. Iniciamos
com um exemplo de matriz  U triangular superior.
3 2 1 0
 0 1 2 3 
U = 0 0 −2 1 


0 0 0 4
Observe que de forma análoga ao caso triangular inferior, U é triangular
superior significa que uij = 0 : j < i. Feita esta observação, vamos deduzir a

7
expressão do termo xi que nos permite construir o algoritmo de resolução. Para
tanto, considere uma linha i fixa do sistema:

n
X
uij xj = yi i = n, n − 1, . . . , 1
j=1
i−1
X n
X
uij xj + uij xj = yi i = n, n − 1, . . . , 1
j=1 j=i
n
X
uii xi + uij xj = yi i = n, n − 1, . . . , 1
j=i+1

Desta forma, a expressão que permite deduzir o algoritmo é dada por:

Pn
yk − j=k+1 ukj xj
xk =
ukk
Cabe destacar que o algoritmo para resolvermos U x = y opera sobre o
sistema linear na ordem inversa das linhas, isto é, primeiro na linha n, depois
na linha n − 1 e assim por diante até trabalhar a linha de ı́ndice 1. Isso porque
para calcular a grandeza xk é necessário dispor das incógnias xn , xn−1 , xk+1 já
calculadas. Por esta razão, o algoritmo é chamado de Algoritmo de Substituições
Retroativas e é apresentado na Figura 3.

function [x] = SubsRetroativas(U,y,n)

for i=n:-1:1
soma = 0.0;
for k = i+1:n
soma = soma + U(i,k)*x(k)
end
if (U(i,i) <> 0)
x(i) = (y(i) - soma)/U(i,i);
else
printf(’Matriz U e singular \n’);
end
end
endfunction

Figura 3: Algoritmo de Substituições Retroativas.

O exemplo abaixo ilustra o uso do algoritmo.

Exemplo 3 --> U,y’

U =

3. 2. 1. 0.
0. 1. 2. 3.
0. 0. -2. 1.

8
 0. 0. 0. 4.
ans =

-10. 10. 1. 12.

--> x = SubsRetroativas(U,y,size(y,1));

--> x’
ans =

-3. -1. 1. 3.

2.2 Resolvendo sistemas lineares a partir de sistemas tri-

angulares
Nesta seção, vamos ilustrar como podemos usar os algoritmos de Substitições
Sucessivas (Figura 1) e Retroativas (Figura 3) para resolver, em duas etapas,
um sistema linear cuja matriz de coeficientes tenha sido fatorada.
Como ainda não apresentamos como produzir a fatoração P A = LU , por
hora, vamos empregar o algoritmo para fatoração P A = LU disponı́vel no
Scilab para obtermos os fatores necessários. Para o exemplo que segue, em-
pregamos a função ResolveParTriangulares(A,b), descrita na Figura 4.

function [x,L,U,P] = ResolveParTriangulares(A,b)

[L,U,P] = lu(A)
[m,n] = size(A)
[y] = SubsSucessivas(L,P*b,n)
[x] = SubsRetroativas(U,y,n)
endfunction

Figura 4: Resolução de um sistema linear Ax = b por meio de dois sistemas

lineares triangulares.

Exemplo 4 Este exemplo ilustra a chamada da funçõ lu do scilab. Considere

então a matriz A e seus fatores.
A =
6. 4. 2. 0.
-3. 0. 3. 6.
9. 7. 7. 2.
12. 9. 12. 12.
b’ =
2. 3. 2. 9.
-->[x,L,U,P] = ResolveParTriangulares(A,b)
x =
1.480D-16
1.0000000
-1.0000000
1.0000000

9
 L =
1. 0. 0. 0.
-0.25 1. 0. 0.
0.5 -0.2222222 1. 0.
0.75 0.1111111 1. 1.
U =
12. 9. 12. 12.
0. 2.25 6. 9.
0. 0. -2.6666667 -4.
0. 0. 0. -4.
P =
0. 0. 0. 1.
0. 1. 0. 0.
1. 0. 0. 0.
0. 0. 1. 0.

Agora vamos usar os fatores para resolver o sistema linear Ax = b, onde

b = (2, 3, 2, 9)T .

3 Fatoração A = LU e P A = LU
3.1 Eliminação de Gauss e Fatoração A = LU
Nesta seção, vamos recordar o método de Eliminação de Gauss, assumindo que
não seja necessário efetuar trocas de linhas do sistema linear. Vamos mostrar
que a Eliminação de Gauss produz os fatores desejados U e L que desejamos.
Por razões didáticas, de inı́cio não faremos uso de trocas de linhas do sistema
linear. Esta é uma hipótese não realista, é adotada aqui apenas para facilitar a
exposição inicial. Com isso, faremos a fatoração A = LU de A e, posteriormente,
ao permitirmos a troca de linhas de A, faremos a fatoração P A = LU .
A ideia da Eliminação de Gauss é transformar o sistema linear Ax = b em
outro sistema linear, U x = y, equivalente ao primeiro. Dois sistemas lineares
equivalentes são indicados como

Ax = b ∼ U x = y,

que significa que toda solução de Ax = b também é solução de U x = y e vice-

versa.
Para transformar Ax = b no equivalente U x = y podemos usar as seguin-
tes operações linha elementares, isto é, operações realizadas sobre as linhas do
sistema que não alteram seu conjunto de soluções:

(T1) Troca da ordem de duas linhas do sistema linear. Como mencionamos,

nesta seção, vamos assumir que não será necessário aplicar T1.
(T2) Multiplicação de uma linha por uma constante não nula.

(T3) Substituição de uma linha do sistema pela soma da própria linha mais um
múltiplo de outra linha do sistema linear.

10
 Sendo i a linha que será substituı́da e j a linha que será multiplicada por
m, temos que
Xn
aik xk = bi
k=1

é substituı́da por
n
X
(aik + majk )xj = bi + mbj .
k=1

Vamos recordar a Eliminação de Gauss por meio de um exemplo. O resultado

será o sistema U x = y, onde este vetor y é exatamente a solução do sistema
linear Ly = b que obterı́amos, caso já dispuséssemos da fatoração A = LU da
matriz. Todas as operações que fizermos sobre as linhas de A, replicaremos nas
linhas de b. Isso é facultativo, pois podemos obter o resultado destas opearações
sobre b, posteriormente, obtendo o y que resolve Ly = b.
Ao longo da aplicação da Eliminação de Gauss, adotaremos a notação de
representar o sistema linear sendo transformado por [Aj |bj ] ∈ Rn×(n+1) , onde
j indica o ı́ndice da operação de pivoteamento completa realizada até aquele
momento do algoritmo (não confundir com a potência j da matriz - não é
o caso aqui). No inı́cio do procedimento, antes de qualquer operação, temos
A0 = A, b0 = b, de forma que [A0 |b0 ] nada mais é do que matriz A expandida
em uma coluna por b. Ao longo da primeira iteração, j = 1, operamos sobre
[A0 |b0 ] e ao final da operação completa obtemos então [A1 |b1 ].
Cada iteração do algoritmo será indexada por uma op j, operação de pivotea-
mento j, que compreenderá um conjunto de operações T2 e T3, para transformar
uma coluna do sistema linear em uma coluna de um sistema triangular superior.
Ao final da j−ésima op completa, temos o sistema equivalente [Aj |bj ]. Faremos
uma op para cada um das colunas de A, exceto a última, visando transformar a
matriz A0 em uma triangular superior: ao final, An−1 será triangular superior.
Ou seja, indexaremos j de 1 até n − 1, inclusive. A ideia de indexação se repete
para as demais iterações.
Em cada operação de pivoteamento faremos uso de um multiplicador mij que
significa: o mutiplicador associado à i−ésima linha de [Aj−1 |bj−1 ], na j−ésima
(op), ou equivalentemente na j−ésima coluna de Aj . Este multiplicador deve
ser escolhido para criar zeros nas linhas abaixo da entrada (Aj−1 )jj . Para um
dado j = 1, . . . , n − 1, calculamos mi,j para i = j + 1, 2, . . . , n.
Para um dado j = 1, . . . , n − 1, mij é dado por:
j−1
aij
mij = − j−1
, i = j + 1, . . . , n. (4)
ajj
Veja que este multiplicador é calculado de forma que a equação abaixo seja
satisfeita:
mij aj−1
jj + aij
j−1
= 0, , i = j + 1, . . . , n. (5)

A linha j de [Aj−1 |bj−1 ] é chamada de pivot. O elemento pivot, aj−1 jj , não

pode ser zero, caso contrário o método falha. Aqui nesta seção, estamos assu-
mindo a hipótese otimista que isso não ocorrerá. Múltiplos desta linha devem
ser somados às linhas de ı́ndice j + 1, j + 2, . . . , n para que a j-ésima coluna de

11
[Aj−1 |bj−1 ] seja transformada em uma coluna de uma triangular superior. Essa
é a ideia da operação T3 que faremos sobre as linhas do sistema.

Exemplo 5 Usando a Eliminação de Gauss (sem pivoteamento parcial, ou

trocade colunas), transformar
  o sistema linear Ax = b em U x = y, para
2 1 1 0 1
 4 3 3 1   3 
A=  8 7 9 5  b =  7 .
  

6 7 9 8 3

op 0 , j = 0 (representação
 do sistema original)

2 1 1 0 1
 4 3 3 1 3 
[A|b] = [A0 |b0 ] = 
 8

7 9 5 7 
6 7 9 8 3
op 1 , j = 1, primeira operação de pivoteamento.
Veja que para j = 1 precisamos multiplicar a linha 1 do sistema por
m2,1 = −2, m3,1 = −4, m4,1 = −3, respectivamente, para criarmos ze-
ros nas posições a02,1 , a03,1 , a04,1 , respectivamente. Empregando estes mul-
tiplicadores, ao finalda primeira operação  de pivoteamento (op) temos o
2 1 1 0 1
 0 1 1 1 1 
sistema: [A1 |b1 ] =  0 3 5 5 3 . A primeira coluna de A é uma
 1

0 4 6 8 0

coluna de uma triangular superior.

 = −3, m4,2 = −4.

op 2 , j = 2. Para esta operação, temos m3,2  Ao final
2 1 1 0 1
 0 1 1 1 1 
da mesma, temos o sistema [A2 |b2 ] =  0 0 2 2 0 . As duas


0 0 2 4 −4
2
primeiras colunas de A são colunas de uma triangular superior.
op 3 , j = 3. Para estaterceira operação, temos
 m4,3 = −1. Ao fina da op 3,
2 1 1 0 1
 0 1 1 1 1 
temos: [A3 |b3 ] =   0 0 2 2 0 .


0 0 0 2 −4
Veja que ao final de n − 1 = 3 operações, obtivemos o sistema linear
triangular superior acima. Para resolver o sistema linear, resolvemos por
 T
substituições sucessivas e encontramos x = −1 1 2 −2

3.2 Reinterpretando a Eliminação de Gauss como um con-

junto de transformações lineares
Cada operação T3, que representa a multiplicação da linha pivot j por mij , a
soma do resultado com a linha i e sua substituição pela soma, pode ser repre-
sentado por um produto de uma matriz Mj por [Aj−1 , bj−1 ]. Atenção aqui: Mj
é a matriz de multiplicadores e não uma coluna desta matriz.

12
 A matriz Mj , chamada de matriz de multiplicadores na j−ésima op, é uma
matriz que difere da matriz identidade apenas pela sua j−ésima coluna, nos
elementos das linhas i : i > j, isto é (abaixo da diagonal principal), que recebem
os multiplicadores mij calculados através de (4). A tı́tulo de ilustração, veja a
j−1
forma da matriz de multiplicadores  na primeira op, para uma matriz A , com
1 0 0 0
 m21 1 0 0 
4 linhas M1 =  m31 0 1 0 


m41 0 0 1
No caso do exemplo que ilustramos, a instanciação da matriz M1 corresponde
à matriz  
1 0 0 0
 −2 1 0 0 
M1 =  −4 0 1 0 .


−3 0 0 1
Observe que a matriz Mj é uma triangular inferior, com diagonal unitária.
Seu determiante corresponde ao produto dos elementos em sua diagoanal, de
forma que det(Mj ) = 1 e a matriz admite inversa.

Exercı́cio 3.1 Verifique que a matriz inversa de Mj , Mj−1 é uma matriz iden-
tifidade, exceto pela j− ésima coluna, que recebe o simétrico das entradas de
Mj , nas linhas
 abaixo da diagonal principal: 
1 0 0 0 0 0 0
 0 1 0 0 0 0 0 
 
 0 0 ...
 
0 0 0 0 
Mj−1 = 
 
 0 0 0 1 0 0 0 

 0 0 0 −mj+1,j 1 0 0 
 
 .. .. 
 0 0 0 . 0 . 0 
0 0 0 −mn,j 0 0 1

Definida a matriz Mj , veja que cada op pode ser represetada da seguinte

forma:

M j [Aj−1 , bj−1 ] = [Aj−1 , bj−1 ]. (6)

No processo de Eliminação de Gauss que ilustramos anteriormente, não fize-
mos isso explicitamente. Nossa ideia é usar a representação (6), para obtermos
os fatores.
Veja o processo completo associado à Eliminação de Gauss que ilustramos.

Exemplo 6 O objetivo deste exemplo é ilustrar, passo a passo, as três trans-

formações lineares que triangularizam a matriz A do Exemplo 5.
    
1 0 0 0 1 0 0 0 1 0 0 0 2 1 1 0 1
 0 1 0 0 
 0 1 0 0   m2,1 1 0 0   4 3 3 1 3 
  
 =
 0 0 1 0   0 m2,3 1 0   m3,1 0 1 0   8 7 9 5 7 
0 0 m4,3 1 0 m2,4 0 1 m4,1 0 0 1 6 7 9 8 3

13
     
1 0 0 0 1 0 0 0 1 0 0 0 2 1 1 0 1
 0 1 0 0   0 1 0 0   −2 1 0 0   4 3 3 1 3 
    =
 0 0 1 0  0 −3 1 0   −4 0 1 0  8 7 9 5 7 
0 0 −1 1 0 −4 0 1 −3 0 0 1 6 7 9 8 3

 
2 1 1 0 1
 0 1 1 1 1 
 .
 0 0 2 2 0 
0 0 0 2 −4

Vamos agora formalizar o processo de Eliminação de Gauss, de forma que

obtenhamos os fatores L, U desejados na fatoração. Veja que realizamos n − 1
operações, uma por coluna de A, e que por construção obtivemos:

Mn−1 Mn−2 . . . M2 M1 [A|b] = [U |y]

⇒
Mn−1 Mn−2 . . . M2 M1 A = U
Mn−1 Mn−2 . . . M2 M1 b = y

Veja o processo de pré-multiplicar A pelas matrizes M1 , M2 , . . . , Mn−1 , nesta

ordem, gerou uma matriz U . Em linguagem de ALC, triangularizamos A. O
mesmo conjunto de transformações lineares aplicados em b gerou o vetor y que
resolve o sistema (1) (até aqui a matriz P é a identidade, pois não trocamos
linhas de A).
Vamos agora considerar a matriz M = Mn−1 Mn−2 . . . M2 M1 . Observe que
M corresponde ao produto de n − 1 matrizes de multiplicadores, cada qual com
determinante igual a 1. Portanto det(M ) = 1n−1 = 1 e M admite inversa. Sua
inversa é dada por
L = M −1 = M1−1 M2−1 . . . Mn−2
−1 −1
Mn−1
A matriz L é uma triangular inferior, com diagonal unitária, que difere
da matriz identidade, pois as entradas abaixo da diagonal principal recebem o
simétrico dos multiplicadores, isto é, Lij = −mij .

Exemplo 7 Verificar que M −1 = L é uma triangular inferior com diagonal

unitária, e suas entradas abaixo da diagonal são os simétricos dos multiplicado-
res dados por (4).
 
1 0 0 0 0
 −m21 1 0 0 0 
 
 . . 
L =  −m31 −m32
 . 0 0 
 (7)
 .. .. . .

 . . . 1 0 
−mn,1 −mn,2 · · · −mn,n−1 1

A verificação se dá analisando os produtos das matrizes que definem L, dos

fatores mais à direita para os mais à esquerda. Veja que o resultado do primeiro
−1 −1
produto Mn−2 Mn−1 :

14
 −1 −1
Mn−2 Mn−1 =
  
1 0 0 0 0 1 0 0 0 0
 ..  0 1 0 0 0 
 0 . 0 0 0   
= 0
  .. 
 0 1 0 0 

 0 0 . 0 0 

 0 0 −mn−1,n−2 1 0   0 0 0 1 0 
0 0 −mn,n−2 0 1 0 0 0 −mn,n−1 1
 
1 0 0 0 0
 .. 
 0 . 0 0 0 
= 0
 
 0 1 0 0 

 0 0 −mn−1,n−2 1 0 
0 0 −mn,n−2 −mn,n−1 1
−1 −1
Observe que a matriz resultante do produto Mn−2 Mn−1 tem a seguinte forma:
• Suas primeiras n − 3 colunas são colunas de uma identidade de ordem n.
A mesma observação se aplica para a última coluna, que também é uma
coluna de uma identidade.

• As demais colunas, de ı́ndices n − 2 a n − 1, são as colunas de ı́ndices

−1 −1
n − 2 e n − 1, respectivamente de Mn−2 e de Mn−1 .
−1 −1
O padrão a ser observado é que a matriz correspondente ao produto Mn−k · · · Mn−1
(para algum k ≥ 2) possui suas primeiras n − k − 1 primeiras colunas, assim
como a última, de ı́ndice n, como colunas da identidade. Além disso, as colu-
nas de ı́ndice n − k até n − 1 são as colunas de ı́ndices (n − k), · · · , (n − 1) de
−1 −1
Mn−k , . . . , Mn−1 , respectivamente.
Assim, ao incorporamos o próximo fator no cálculo de L, isto é, ao calcu-
−1 −1 −1
larmos Mn−k−1 (Mn−k · · · Mn−1 ) preservamos as primeiras n − k − 1 colunas de
−1
Mn−k−1 , incluindo sua coluna de ı́ndice n − k − 1, que é a única de suas colunas
que difere da identidade e que armazena os simétricos dos multiplicadores da op
−1
n − k − 1. Todas as demais colunas de Mn−k−1 são colunas da identidade. Por-
−1 −1 −1
tanto, o produto Mn−k−1 (Mn−k · · · Mn−1 ) para as colunas de ı́ndice n − k em
diante será uma combinação linear das colunas da identidade por pesos que vêm
−1 −1 −1 −1 −1
das colunas de (Mn−k · · · Mn−1 ). Já as colunas de Mn−k−1 (Mn−k · · · Mn−1 ) de
ı́ndice n−k−1 ou menor serão colunas da identidade. Esta é a invariante do pro-
cesso, que se repete até incorporarmos o fator M1−1 ao produto M2−1 . . . Mn−1 −1
,
para k = n − 2.
Para ilustrar a incorporação de mais um fator aos já avaliados, considere o

15
 −1 −1 −1
resultado de Mn−3 Mn−2 Mn−1 :
−1 −1 −1
Mn−3 Mn−2 Mn−1 =
1 0 0 0 0 0 1 0 0 0 0 0
  
 ..  .. 

 0 . 0 0 0 0 
 0 . 0 0 0 0 

=
 0 0 1 0 0 0 
 0 0 1 0 0 0 


 0 0 −mn−2,n−3 1 0 0 
 0 0 0 1 0 0 

 0 0 −mn−1,n−3 0 1 0  0 0 0 −mn−1,n−2 1 0 
0 0 −mn,n−3 0 0 1 0 0 0 −mn,n−2 −mn,n−1 1
1 0 0 0 0 0
 
 .. 

 0 . 0 0 0 0 

=
 0 0 1 0 0 0 


 0 0 −mn−2,n−3 1 0 0 

 0 0 −mn−1,n−3 −mn−1,n−2 1 0 
0 0 −mn,n−3 −mn,n−2 −mn,n−1 1

Em resumo, desde que armazenados os multiplicadores empregados no pro-

cesso, a Eliminação de Gauss produz os fatores L, U de A, onde L é triangular
inferior com diagonal unitária e U é triangular superior. Veja:

Mn−1 Mn−2 . . . M2 M1 [A|b] = [U |y]

A = M1−1 M2−1 . . . Mn−2
−1 −1
Mn−1 U
M1−1 M2−1 . . . Mn−2
−1 −1
Mn−1 =L
A = LU

3.3 Visão por colunas da Fatoração A = LU

Uma maneira bastante conveniente de se formalizar as operações da Eliminação
de Gauss é verificar que, a cada op, subtraı́mos de A uma matriz de rank 1, e
depois operamos sobre a diferença, repetindo o processo, até a última op. Esta
visão alternativa é chamada de visão coluna, ou visão de soma de matrizes de
posto 1, para a fatoração.
Para deduzirmos esta visão alternativa, vamos partir da fatoração A = LU
que resulta da Eliminação de Gauss, escrevendo-a como uma soma de n =
posto(A) matrizes de posto 1.

A = LU
uT1
 
 
| | ... |  uT2 
A =  L1 L2 ... Ln   
 ... 
| | ... |
uTn
n
X
= Lj uTj ,
j=1

onde Lj corresponde à j−ésima coluna de L e uTj é a j−ésima linha de U . O

último termo da soma, Ln uTn não foi explicitamente calculado na Eliminação
de Gauss, pois não era necessário (já que a última coluna de An−1 sempre é
uma coluna de uma triangular superior de ordem n), sendo trivialmente dado

16
por en uTn (recorde-se de nossa notação que utiliza ei para representar um vetor
n−dimensional de zeros, exceto pela i−ésima entrada que é 1). Veja que para
todo ı́ndice j = 1, . . . , n, as linhas uTj são as linhas pivotais: a linha j da matriz
Aj−1 obtida ao longo da Eliminação de Gauss. Já as colunas de Lj satisfazem:
lij = 0 para i < j, ljj = 1, lij = −mij .
Naturalmente, quando aplicarmos a Eliminação, não temos todos os termos
Lj uTj para j = 1, . . . , n. Estes termos são descobertos ao longo do processo.
Porém, veja que para um determinado ı́ndice j de op, por exemplo, j = 1,
temos

n
X
A= Lk uTk
k=1
n
X
= L1 uT1 + Lk uTk
k=2
n
X
A − L1 uT1 = Lk uTk
k=2
= A2
Pn
A matriz A2 = k=2 Lk uTk é uma soma de n − 2 matrizes de posto 1. Nova-
mente salientamos que, só ao final da Eliminação de Gauss é que dispomos dos
demais termos L2 uT2 , . . . , Ln uTn na fatoração. Mas podemos aplicar a mesma
ideia agora à matriz A2. Veja que a matriz A2 corresponde às últimas n − 1
colunas e linhas da matriz A1 que obtivemos ao final da primeira op na Eli-
minação de Gauss. Então, a linha uT2 é a linha 2 de A1 e a coluna L2 é obtida
calculando-se os multiplicadores pertinentes à segunda op. Incorporando mais
uma op ao processo temos:

n
X
A − L1 uT1 = Lk uTk
k=2
= A2
Xn
A − L1 uT1 − L2 uT2 = Lk uTk
k=3
= A3
Pj−1
Então, a matriz Aj para j = 1, . . . , n é simplesmente A − k=1 Lk uTk , onde
A1 = A. Repetindo o processo por n − 1 ops, temos que
n−1
X
A− Lj uTj = en uTn ,
j=1

Pn−1
o que nos permite escrever A = j=1 Li uTi + en uTn , sendo a fatoração final.
Exemplo 8 Vamos ilustrar a visão de colunas da Fatoração, por meio do exem-
plo da seção anterior, interpretando agora cada op como a subtração de uma

17
matriz de posto 1 que é dada pelo produto externo de uma coluna de L por uma
linha de U (a linha pivotal de Aj−1 ).

1. primeira op
   
1 vez a linha pivot 1 1 vez a linha pivot 1
 l21 vez a linha pivot 1  l21 vez a linha pivot 1 
• A=

+A2 =  +
 l31 vez a linha pivot 1   l31 vez a linha pivot 1 
l41 vez a linha pivot 1 l41 vez a linha pivot 1
 
0 0 0 0
 0 × × × 
 
 0 × × × 
0 × × ×

• Recorde-se dos valores que foram calculados para os multiplicadores

que permitem escrever a primeira coluna de L com as entrada l21 =
a21 /a11 , l31 = a31 /a11 , l41 = a41 /a11 . Recorde-se também que as
entradas l21 , l31 , l41 são os simétricos de m21 , m31 , m41 .
• 
Então temos      
2 1 1 0 1 0 0 0 0 2 1 1 0
 4 3 3 1   2    0 1 1 1   4 3 3 1 
 =  2 1 1 0 +  =
 8 7 9 5   4   0 3 5 5  8 7 9 5 
6 7 9 8 3 0 4 6 8 6 7 9 8
   
2 1 1 0 0 0 0 0
 4 2 2 0   0 1 1 1 
 + 
 8 4 4 0   0 3 5 5 
6 3 3 0 0 4 6 8
2. segunda op:
• A2 = L2 uT2 + A3.
     
0 0 0 0 0 0 0 0 0
 0 1 1 1   1   0 0 0 0 
• A2 = 

=  0 1 1 1 + 
 0 3 5 5   3   0 0 2 2 
0 4 6 8 4 0 0 2 4
     
0 0 0 0 0 0 0 0 0 0 0 0
 0 1 1 1   0 1 1 1   0 0 0 0 
• A2 = 
 0
= + 
3 5 5   0 3 3 3   0 0 2 2 
0 4 6 8 0 4 4 4 0 0 2 4
3. terceira op temos:

• A3 = L3 uT3 + A4.
     
0 0 0 0 0 0 0 0 0
 0 0 0 0   0   0 0 0 0 
• A3 = 

=  0 0 2 2 + 
 0 0 2 2   1   0 0 0 0 
0 0 2 4 1 0 0 0 2
4. o quarto termo, referente à uma quarta op não necessária na visão por
linhas, mas necessária na visão colunas, corresponde a:

18
 
0
 0 
• A4 = L4 uT4 = e4 uT4 = 

 0  0 0 0 2


3.4 Complexidade computacional de A = LU

Na Figura 5, apresentamos o Algoritmo que produz a fatoração A = LU , à
partir da Eliminação de Gauss, sem ainda incorporar a troca de linhas. Vamos
discutir sua complexidade computacional e mostrar que é um algoritmo na classe
O(n3 ). O algoritmo usa a notação B(p : k, l : m) para representar a submatriz
de B que contém as linhas de ı́ndices p até k e colunas de ı́ndices l até m de B.

function [U,L] = EliminacaoGauss(A,n)

U = A
L = eye(n,n)
for j=1:n-1
for i = j+1:n
L(i,j) = U(i,j)/U(j,j)
U(i,j:n) = U(i,j:n) - L(i,j)*U(j,j:n)
end
end
endfunction

Figura 5: Algoritmo para Fatoração A = LU .

A instrução mais relevante para a complexidade computacional é

U(i,j:n) = U(i,j:n) - L(i,j)*U(j,j:n)

Observe que o termo

L(i,j)*U(j,j:n)
na instrução em estudo corresponde ao produto de um escalar por um vetor n−j
dimensional. Portanto, requer n − j operações aritméticas de ponto flutuante.
Vamos analisar apenas P a complexidadePnadicionada por esta instrução, que é
n−1 Pn
dada então pela soma j=1 i=j+1 k=j+1 1. Para o desenvolvimento que
segue, recorde-se do valor das somas
n
X n(n + 1)
i=
i=1
2

e também
n
X n
i2 = (n + 1)(2n + 1).
i=1
6
Então, a função que determina o número de vezes que a instrução

U(i,j:n) = U(i,j:n) - L(i,j)*U(j,j:n)

19
 será executada pode ser obtida desenvolvendo-se

n−1
X n
X n
X
1 =
j=1 i=j+1 k=j+1
n−1
X n
X
(n − j) =
j=1 i=j+1
n−1
X
(n − j)2 =
j=1
n−1
X n−1
X n−1
X
n2 − 2n j+ j2 =
j=1 j=1 j=1
n−1
X
j2 =
j=1
n−1 2n3 −3n2 +n
n(2n − 1) = 6
6
3 2
Veja que a complexidade adicionada pela instrução é 2n −3n 3
+n
, pois cada
vez que a instrução for executada, serão realizadas uma subtração e uma soma.
Portanto, a fatoração A = LU custa O(n3 ) operações aritméticas, e a constante
do termo cúbico na função de complexidade é 23 .

3.5 Introduzindo o pivoteamento de colunas: P A = LU

3.5.1 Problemas ao não se trocar as linhas da matriz
A Eliminação de Gauss na forma como apresentamos até aqui, sem incorporar
a troca de linhas de A, não é prática e não funciona para a quase totalidade dos
casos de interesse. Foi apresentada apenas por razões didáticas, visando ilustrar
que, se forem armazenados os multiplicadores e se não houver divisão por zero
em (4), a Eliminação produz os fatores L e U de A.
O problema da Eliminação de Gauss, no entanto, não se resume ao caso em
que aj−1
jj em (4) é identicamente nulo, caso em que j−ésima iteração não seria
definida no algoritmo da Figura (5) que apresentamos. Quando o denominador
em (4) não é zero, mas é muito pequeno comparado ao numerador, o módulo
do multiplicador tende a ser muito grande. Esse multiplicador será utilizado
nas transformações lineares (5) e, desta forma, erros numéricos muito represen-
tativos devem ser observados. Recorde-se das duas fontes principais de erros
numéricos:
• Subtração de quantidades muito próximas.
• Soma de quantidades muito dı́spares.
10−17 1

Exemplo 9 Para este exemplo, seja A = a matriz a ser fato-
1 1
rada. Veja que de acordo com (4) e (5), 1017 vezes a primeira linha é sub-
traı́do da segunda linha. Assumindo que utilizemos aritmética de precisão infi-
nita, o que não é o caso com o uso de computadores digitais, a Eliminação de

20
Gauss produziria os seguintes
 −17 fatores exatos(sem erros numéricos) L, U abaixo.
1 0 10 1
L= ,U=
1017 1 0 1 − 1017
Considere agora uma condição realista, em que empregamos aritmética de
precisão finita e a precisão da máquina é ϵ ≈ 10−16 . A grandeza 1 − 1017 na
entrada de L não seá representada de forma exata. Ao invés disso, obteremos
o resultado −1017  . Desta forma, os fatores obtidos com ϵ ≈ 10−16 são L̃, Ũ
−17
 
1 0 10 1
dados por L̃ = , Ũ = , cujo produto é L̃Ũ =
1017 1 0 −1017
 −17 
10 1
, uma matriz substancialmente distinta de A. Se agora desejarmos
1 0
resolver o sistema linear Ax = b para b = (1, 0)T via L̃Ũ x = b, obteremos
x̃ = (0, 1)T . Porém, a solução verdadeira do sistema linear é x = (−1, 1)T . Ou
seja, a solução numérica e a solução verdadeira são muito distantes. Verifique
o resultado deste experimento numérico no seguinte código scilab.
-->Ltilde = [1 0;1E17 1]
Ltilde =
1. 0.
1.000D+17 1.
-->Utilde = [1E-17 1;0 1-1E17]
Utilde =
1.000D-17 1.
0. -1.000D+17
-->[y] = SubsSucessivas(Ltilde,b,2)
y =
1.
-1.000D+17
-->[x] = SubsRetroativas(Utilde,y,2)
x =
0.
1.
Entretanto, observe que implementarmos o pivoteamento parcial, teremos o re-
sultado correto. Veja o resultado com o uso da função lu do scilab.
-->[x,L,U,P] = ResolveParTriangulares(A,b)
x =
-1.
1.
L =
1. 0.
1.000D-17 1.
U =
1. 1.
0. 1.
P =
0. 1.
1. 0.
A resposta do procedimento ResolveParTriangulares(A,b), ilustrado na Fi-
gura 4, é composta pela solução x, a matriz de permutação P e os fatores L, U ,

21
de forma que P A = LU . A matriz P mostra que, alterando a ordem das li-
nhas da matriz A, o problema foi resolvido. Observe que, diante da troca de
linhas, os fatores L, U foram calculados sem que os erros numéricos inerentes
à computação digital produzissem uma resposta muito diferente da verdadeira.

Pelas razões discutidas acima e ilustradas no exemplo, verificamos que os

multiplicadores na Eliminação de Gauss não devem ser muito grandes. Por
meio da operação de pivotemento de linhas, isto é, da troca de linhas de A,
garantimos que os multiplicadores empregados tenham módulo não superior a
1. Como produzir esta fatoração é o assunto da próxima seção.

3.5.2 Fatoração P A = LU
O instrumento que usaremos para representar a troca de linhas de A visando
reduzir os erros numéricos e divisão por zero é o de matriz de permutação.
Dada uma permutação π = (π1 , π2 , . . . , πn ) dos inteiros {1, . . . , n}, uma matriz
P é uma permutação (de linhas, associado a π) da matriz identidade In se a
i-ésima linha de P for a πi -ésima linha de In . Duas propriedades importantes
das matrizes de permutação devem ser recordadas aqui. Se P é matriz de
permutação, sua inversa é sua transposta, isto é: P −1 = P T . Uma matriz de
permutação é um caso particular de uma ortogonal (ou unitária). A segunda
propriedade é que o determinante de uma matriz de permutação é (−1)p , onde
p é o número de trocas de linhas necessárias realizadas em P , para que a matriz
resultante destas trocas seja a identidade de mesma ordem.
Defina  como uma matriz que contém as mesmas linhas de A, apenas
apresentadas em ordem diferente. Assuma que as linhas em que são apresentadas
em  seja dada por uma permutação π dos ı́ndices das linhas de A. Então
existe uma matriz de permutação P tal que  = P A. Se A é não singular (logo,
det(A) ̸= 0), então existe uma matriz de permutação P tal que podemos aplicar
o Método de Eliminação de Gauss à matriz P A, sem que ocorra divisão por
zero, no cálculo dos multiplicadores. Logo P A = LU . Assim, a permutação de
linhas resolve o primeiro problema que identificamos, que é a divisão por zero.
Resta-nos eleger algum bom critério para permutar as linhas de A, obtendo P e
P A. Essa matriz P será descoberta ao longo do processo, ao longo da Eliminação
de Gauss.
O critério para definir uma boa P é o seguinte. Desejamos uma P tal que
a matriz L obtida ao se fatorar P A tenha entradas cujos módulos sejam no
máximo iguais a 1. Para garantir esta propriedade, em cada op de ı́ndice j,
a linha pivotal não será necessariamente a linha j de Aj−1 . Vamos comparar
as entradas na coluna j, nas linhas j, j + 1, . . . , n de Aj−1 , e eleger como linha
pivotal p aquela que contiver o maior elemento em módulo naquela coluna. O
pivoteamento é chamado de parcial pois envolve a comparação dos módulos
apenas nas linhas e não na matriz toda, sem que haja troca de colunas de
A também. O pivoteamento total é uma alternativa mais cara: escolhemos
o elemento de maior módulo da submatriz quadrada de Aj−1 que envolve as
colunas e linhas de j até n. Isso implicaria em trocar a ordem de colunas e linhas
e pesquisar o maior elemento dentre O((n − j + 1)2 ) alternativas, elevando o
custo total. Também seriam necessárias duas matrizes de permutação, uma para
troca de linhas, P , e outra para troca de colunas, P̂ , de forma que P AP̂ = LU .
A permutação total não será empregada aqui.

22
 Resumindo então o que fazemos na permutação parcial, com troca de linhas:
na j−ésima op, escolhemos como elemento pivô o valor

aj−1 j−1
p,j = max{|ak,j | : k = j, . . . , n}

e a linha pivotal como a linha

j−1
p = arg max{|ak,j | : k = j, . . . , n}.

Caso p ̸= j, trocamos as linhas p e j de Aj−1 antes de fazer as operações. No

exemplo e no algoritmo que apresentaremos na sequência, não armazenaremos
explicitamente a matriz P , mas sim um vetor auxiliar, pivot, que guarda a
ordem das linhas pivotais. O valor inteiro pivot(i) indica o ı́ndice da linha
do sistema original representado na linha i. Este vetor será inicializado como
pivot(i) = i, i = 1, . . . , n, assumindo que não haverá troca de linhas. Sempre
que alguma troca ocorrer, trocamos o conteúdo armazenado em pivot(j) por
pivot(p), e vice versa. Vamos ilustrar o processo por meio de um exemplo.

Exemplo 10 Vamos resolver o sistema linear apresentado no Exemplo 5, usando

a fatoração P A = LU . Vamos também usar a fatoração para o cálculo do de-
terminante. Neste exemplo, indicamos as linhas pivotais pelo elemento em ver-
melho de maior módulo nas linhas que competem para serem as linhas pivotais.
A ordem das linhas do sistema original é indicada à direita de [Aj |bj ]. Esta
ordem deve refletir as informações armazenadas no vetor pivot. Agora, com o
pivoteamento parcial, podemos ter multiplicadores do tipo mii (com ı́ndices de
coluna e linha iguais), pois a linha i pode não ter sido a linha pivotal na op i.
Diante desta abordagem, o primeiro ı́ndice i associado ao multiplicador mij não
faz referência à posição fı́sica i da linha, mas sim qual linha i do sistema original
é representada na posição considerada para o cálculo dos multiplicadores.

1. Inicialização: pivot = (1, 2, 3, 4)T .

 
2 1 1 0 1 (E1)
 4 3 3 1 3 (E2) 
2. (op1) j = 1, [A|b] =   
8 7 9 5 7 (E3) 
6 7 9 8 3 (E4)
• ap1 = 8, p = 3.
• Atualizamos “pivot(1) = pivot(3), pivot(3) = pivot(1)”, trocamos o
conteúdo das linhas 1/3.
• Pre-multiplicamos por P1 , que difere de I nas linhas 1 e 3 apenas.
• m11 = − 14 , m21 = − 12 , m41 = − 34
 
8 7 9 5 7 (E3)
 0 −1 −3 −3 −1
2 2 2 2 (E2′ ) 
M1 P1 [A|b] = 
(E1′ ) 

 0 −3 −5 −5 −3
4 4 4 4
0 7
4
9
4
17
4 − 94 (E4′ )
 
8 7 9 5 7 (E3)
 0 −1 −3 −3 −1
2 2 2 2 (E2′ ) 
3. (op2) j = 2, M1 P1 [A|b] = 
(E1′ ) 

 0 −3 −5 −5 −3
4 4 4 4
0 7
4
9
4
17
4 − 94 (E4′ )

23
 • ap2 = 74 , p = 4.
• Trocamos o conteúdo da linha j = 2 pela linha p = 4, que corresponde
a fazer “pivot(2) = pivot(4), pivot(4) = pivot(2)”.
• Pré-multplicamos por P2 , que difere de I nas linhas 2 e 4 apenas.
• m22 = 72 , m12 = 3
7
 
8 7 9 5 7 (E3)
 0 7
4
9
4
17
4 − 49 (E4′ ) 
M2 P2 M1 P1 [A|b] = 
(E1′′ ) 

 0 0 − 27 4
7 − 12
7
0 0 − 67 − 27 − 78 (E2′′ )
 
8 7 9 5 7 (E3)
 0 7
4
9
4
17
4 − 94 (E4′ ) 
4. (op3) j = 3, M2 P2 M1 P1 [A|b] = 
(E1′′ ) 

 0 0 − 72 4
7 − 12
7
0 0 − 67 − 27 − 87 (E2′′ )

• ap3 = − 67 , p = 4.
• Trocamos o conteúdo da linha j = 3 pela linha p = 4, que corresponde
a fazer “pivot(3) = pivot(4), pivot(4) = pivot(3)”.
• Pré-multplicamos por P3 , que difere de I apenas nas linhas 3 e 4.
• m13 = − 13
 
8 7 9 5 7 (E3)
 0 7
4
9
4
17
4 − 94 (E4′ ) 
M3 P3 M2 P2 M1 P1 [A|b] = 
(E2′′ ) 

 0 0 − 76 − 27 − 87
0 0 0 2
3 − 43 (E1′′′ )
Agora, vamos resolver o sistema linear Ax = b. Observe que a solução de
Ly = P b já é disponı́vel, pois operamos sobre [A|b] e não apenas sobre b. Bas-
taria, portanto, resolvermos U x = y. Entretando, vamos explicitar todos os
fatores obtidos com a fatoração e a resolução dos dois sistemas lineares trian-
gulares, recalculando y.
Os fatores já disponı́veis são P (pois dispomos de pivot) e U .
 
0 0 1 0
 0 0 0 1 
• pivot = 3 4 2 1 → P = 
 

 0 1 0 0 
1 0 0 0

 
8 7 9 5
 0 7 9 17
• U =

4 4 4 
 0 0 − 76 − 27 
2
0 0 0 3

Para termos a fatoração completa, o fator que nos resta determinar é L.

Como definir L à partir dos multiplicadores ? Veja que o ı́ndice da linha de A
que gera uma linha pivotal na op j é a linha pivot(j). Portanto, para a linha
pivot(j) de L teremos j −1 multiplicadores, calculados nas ops anteriores. Estes
serão armazenados na linha j de L, na ordem em que foram gerados.

24
 • m11 = − 14 , m21 = − 12 , m41 = − 43 , m22 = 27 , m12 = 37 , m13 = − 31 .

 
1
3
1
• L=
 
4 
1

2 − 27 1 
1
4 − 37 1
3 1

Para resolver o sistema, resolvemos dois sistemas lineares triangulares, já

que Ax = b → P Ax = P b → LU x = P b.
      
1 y1 7 7
 3  9 
1  =   , y =  − 48 
  y2   3 
• Ly = P b :  4
1 2


2 −7 1   y3   3   − 
7
1 3 1
4 − 7 3 1 y 4 1 − 43
      
8 7 9 5 x1 7 −1
 0 7 9 17   x2   − 9
• Ux = y :  , x =  1 
  
4 4 4   4
  x3  =  − 8
6

 0 0 −7 − 72 7
  2 
2
0 0 0 3
x4 − 43 −2

A fatoração também nos permite calcular o determinante da matriz A: det(P A) =

det(LU ) = det(U ). Logo, det(A) = det(U )(−1)p , onde p é o número de tro-
cas de pares de linhas necessárias para transformar P em I. Com isso temos:
det(A) = (8)( 47 )( −6 2 3
7 )( 3 )(−1) = 8.

3.5.3 Visão por colunas da fatoração P A = LU

Para ilustrar como podemos representar a Eliminação com pivoteamento parcial,
através da visão de colunas, vamos escrever A = L̃U onde L̃ não é triangular
inferior como desejamos, mas poderá ser transformada em L, através da per-
mutação de linhas pertinente. Ou seja, vamos escrever uma fatoração A = L̃U ,
onde não vamos nos preocupar com a forma de L̃. Esta matriz apenas deve
representar as transformações lineares que desejamos fazer. As linhas pivotais
que descobrirmos irão abastecer as linhas de U , na ordem em que forem desco-
bertas. Com o processo, vamos descobrir uma matriz de permutação P tal que
P A = (P L̃)U = LU , onde P L̃ é a L que desejamos. Vamos ilustrar o processo
com o exemplo seguinte.

Exemplo 11 • j = 1, primeira op.

Seguindo a mesma estratégia de escolher o pivot de maior módulo, a pri-
meira linha de U é a terceira linha de A. Os (simétricos dos) multiplicado-
res são calculados sem relação a esta linha, obtendo a primeira
 coluna de
0 1 1
L̃. Assim, a primeira op deve produzir o resutlado: A =  1 3 7  →
2 4 8
 
0 1 1
 0 1 3 . Vamos guardar os vetores L̃i sem nos preocupar com a
2 4 8

25
  
  0
forma triangular (inferior) para eles. uT1 = 2 4 8 L̃1 =  0.5 
  1
0 1 1
A = L̃1 uT1 + A2 → A2 =  0 1 3 
0 0 0

• Note que a terceira linha de A2 (e não a primeira como na aplicação do

método sem troca de linhas) é toda composta de zeros. pivot(1) = 3.
• j = 2, segunda op. Há empate para escolha do pivot. Adotamos pivot(2) =
1.
 A segunda operação
 de pivoteamento
 deve gerar o resultado A2 =
0 1 1 0 1 1  
 0 1 3  →  0 0 2 . Então temos: uT2 = 0 1 1 , L̃2 =
0 0 0 0 0 0
   
1 0 0 0
 1 . A2 = L̃2 uT2 + A3 → A3 =  0 0 2 .
0 0 0 0
 
0
• j = 3, terceira op. pivot(3) = 2. uT3 = 0 0 2 , L̃3 =  1 .
 

  0
0 0 0
A3 =  0 0 2 
0 0 0
Veja que neste momento, já podemos escrever a A como soma de matrizes
de rank-1 que acumulamos:
  
0 1 0 2 4 8
• A =  0.5 1 1  0 1 1 
1 0 0 0 0 2

• A = L̃1 uT1 + L̃2 uT2 + L̃3 uT3 = L̃U

A forma acima ainda não é a desejada pois as colunas L̃i não são colunas
de uma triangular inferior com diagonal unitária. Como dispomos de pivot,
calculamos:
 
0 0 1
• pivot = (3, 1, 2) → P =  1 0 0 
0 1 0

• Pré-multiplicando por P a fatoração acima, temos: P A = (P L̃)U = LU

    
2 4 8 1 0 0 2 4 8
• PA =  0 1 1  =  0 1 0  0 1 1 
1 3 7 0.5 1 1 0 0 2

4 Fatoração de Cholesky
A fatoração de Cholesky A = LLT é uma forma particular da fatoração A = LU ,
na qual L = U T e as entradas na diagonal de L são positivas. Isto é, L é

26
triangular inferior e lii > 0 para todo i = 1, . . . , n. A fatoração é possı́vel se e
somente se a matriz A ∈ Sn++ , isto é, se A é simétrica e positiva definida.
A caracterização da positividade de uma matriz simétrica pode ser realizada
de diversas formas. Uma delas consiste em empregar o algoritmo de fatoração de
Cholesky apresentado nesta seção. Se o algoritmo for bem sucedido, chegando
ao final sem efetuar divisão por zero ou radiciação de argumento negativo, o
fator L é produzido e a positividade é caracterizada. Caso contrário, conclui-se
que a fatoração não é possı́vel e que a matriz não é positiva.
Duas são as vantagens da fatoração de Cholesky em relação à Fatoração LU .
A positividade da matriz garante não haver divisão por zero, ou pivot nulo, de
forma que a fatoração é bastante estável numericamente. A segunda é que, em
função de se explorar a simetria de A, o número de operações aritméticas en-
volvidas é aproximadamente a metade das necessárias na Eliminação de Gauss.
De igual forma, há apenas um fator a ser armazenado, L, de forma que a com-
plexidade de memória também é a metade da fatoração LU .
O algoritmo de Fatoração de Cholesky é baseado no seguinte Teorema.

Teorema 4.1 Teorema de Cholesky. Uma matriz simétrica A ∈ Rn×n é po-

n
sitiva definida (A ∈ S++ ) se e somente se possui uma fatoração (chamada
fatoração de Cholesky) da forma:

A = LLT

onde L ∈ Rn×n é uma matriz triangular inferior com diagonal positiva. Esta
fatoração é única.

Uma interpretação do uso do algoritmo para a caracterização da positividade

pode ser dada da seguinte forma. O algoritmo assume que A seja simétrica
positiva definida e usa a definição A = LLT para se obter os fatores envolvidos.
Em uma perspectiva otimista, se tudo der certo com a aplicação do algoritmo,
obtemos os fatores e caracterizamos a positividade. Tudo de uma vez só.
Assim, a ideia é inicialmente assumir que A seja simétrica e positiva definida.
Então, pelo teorema, existe L com a diagonal positiva tal que A = LLT . Pela
definição dos elementos de A, como o produto interno das linhas de L pelas
colunas de LT (ou pelas próprias linhas de L, já que L = LT ) temos que

aij = liT lj : i, j = 1, . . . , n

onde li denota a i-ésima linha de A. A ideia é então percorrer os elementos de

A por colunas, das de menor ı́ndice j para as de maior ı́ndice e, para um ı́ndice
j de coluna fixo, percorrer os elementos das linhas de ı́ndices j = i, i + 1, . . . , n,
calculando as entradas de L. Mais precisamente, para cada par i, j de linha e
coluna de A, usamos a definição de aij para calcular a entrada lij de L. Para
tanto, vamos empregar as definições de ljj e de lji , obtidas por meio da definição
aij = liT lj .
Para definir a expressão analı́tica de ljj , lji , vamos nos recordar que L é
triangular inferior ⇐⇒ lij = 0, j > i. Então:

27
 ajj =ljT lj
X n
2
= ljk
k=1
j
X n
X
2 2
= ljk + ljk
k=1 k=j+1
j
X
2
= ljk
k=1
j−1
X
2 2
=ljj + ljk
k=1

Adotamos a raiz positiva e então obtemos a expressão:

v
u
u j−1
X
ljj = +tajj − 2
ljk para j = 1, . . . , n. (8)
k=1

Veja que o método pode falhar (e com ele a hipótese de positividade) quando
ljj ≤ 0. Se não for esse o caso, para um determinado ı́ndice i > j, temos:

aij =liT lj
X n
= lik ljk
k=1
j
X n
X
= lik ljk + lik ljk
k=1 k=j+1
j
X
= lik ljk
k=1
j−1
X
=lij ljj + lik ljk
k=1

o que nos permite escrever:

Pj−1
aij − k=1 lik ljk
lij = para i = j + 1, . . . , n. (9)
ljj

Assim sendo, a ideia pode ser sistematizada da seguinte forma. Para todo
2
j = 1, . . . , n, calculamos o elemento ljj de acordo com (8). Se ljj > 0, calculamos
os elementos lij , para todos os valores de i = j+1, . . . , n, de acordo com (9). Esta
ideia é sistematizada na implementação em scilab, apresentada na Figura 6.
Algumas implementações de Cholesky, por exemplo aquela disponı́vel no pacote
scilab retorna a matriz triangular superior. A implementação apresentada na
Figura 6 retorna a triangular inferior.

28
function [L] = Cholesky(A)
n = size(A,1)
L = zeros(n,n)
for j = 1:n
soma = A(j,j)
for k = 1:j-1
soma = soma - L(j,k)*L(j,k)
end
if soma > 0.0
L(j,j) = sqrt(soma)
else
print("Matriz nao e SPD \n")
end
for i = j+1:n
soma = A(i,j)
for k = 1:j-1
soma = soma - L(i,k)*L(j,k)
end
L(i,j) = soma / L(j,j)
end
end
endfunction

Figura 6: Algoritmo de Fatoração de Cholesky que explora a definição dos

elementos aij de A para o cálculo dos fatores.

A implementação apresentada na Figura 6 faz uso das expressões (8) e (9),

que foram derivadas através da definição de aij = liT lj . Por esta razão é de-
nominada implementação por linhas ou por produto interno da Fatoração de
Cholesky.
Verifique você mesmo a corretude da fatoração de Cholesky apresentada no
exemplo abaixo.

Exemplo 12 Vamos empregar o algoritmo para verificar a positividade e obter

os fatores L na fatoração de Cholesky de A.
A =
1. -1. 3. -4.
-1. 5. -1. 2.
3. -1. 14. -9.
-4. 2. -9. 22.
-->L = Cholesky(A)
L =
1. 0. 0. 0.
-1. 2. 0. 0.
3. 1. 2. 0.
-4. -1. 2. 1.

29
4.1 Complexidade da Fatoração de Cholesky
A instrução mais custosa do Algoritmo de Cholesky apresentado na Figura 6 é

soma = soma - L(i,k)*L(j,k)

Cada vez que a instrução é executada, duas operações aritméticas de ponto

flutuante são realizadas. Vamos calcular a contribuição desta instrução para a
complexidade total do algoritmo e mostrar que o algoritmo também é O(n3 ),
assim como a fatoração LU .
O número de vezes em que a instrução é chamada é dada pelo somatório
abaixo avaliado.

n
X n
X j−1
X n
X n
X
1= (j − 1)
j=1 i=j+1 k=1 j=1 i=j+1
Xn
= (n − j)(j − 1)
j=1
Xn
= ((n + 1)j − n − j 2 )
j=1

n(n + 1)2 n
= − n2 − (n + 1)(2n + 1)
2 6
n3 2
= + O(n )
6
Pn Pn
Para a dedução acima, fizemos uso de i=1 i = n(n+1)2 e de i=1 i2 = n6 (n+
1)(2n+1). Veja que o total de operações aritméticas incorridas na instrução deve
ser multiplicado por 2. Portanto, a constante do termo cúbico na complexidade
adicionada pela instrução é 31 .
Quando comparada à complexidade da fatoração A = LU , discutida na
Seção 3.4, verificamos que o termo cúbico na função de complexidade de Cho-
lesky é 13 , enquanto que em LU é 23 . Assim sendo, embora sejam assinton-
ticamente equivalentes (ambos são O(n3 )), para um dado valor de n, o custo
computacional da fatoração de Cholesky é aproximadamente a metade do custo
da Fatoração LU .

4.2 Visão por colunas ou outer Cholesky

Da mesma forma como apresentamos uma formulação mais abstrata para a
fatoração LU , na forma de produtos externos, vamos proceder para a Fatoração
de Cholesky.
Veja que se a matriz A a ser fatorada é positiva definida, a11 > 0 deve valer.
Vamos então particionar A em dois blocos, o primeiro deles sendo a11 , uma
matriz positiva definida de dimensão 1 e o bloco de dimensão n − 1, K = A(2 :
n, 2 : n).
Feito o particionamento, escrevemos a matriz A simétrica como o produto
de duas matrizes, triangulares em blocos (ainda não triangulares, apenas trian-

30
gulares em blocos):
wT sT
    
a11 r11 0 r11
A= = (10)
w K s R̂T 0 R̂
Veja que podemos calcular r11 e s facilmente. Na verdade, o conjunto de
operações que apresentamos anteriormente (na versão produto escalar ou linha
da Fatoração de Cholesky) para computar a primeira coluna de L ou (primeira
linha de LT ) corresponde a calcular estas entradas. Veja:

√
r11 = a11
w
s= .
r11
Além disso, sabemos que
R̂T R̂ = K − ssT , (11)
já que o segundo bloco K de A é o produto da segunda linha de blocos do
primeiro fator pela segunda coluna de blocos no segundo fator em (10): K =
R̂T R̂ + ssT . Veja que, ao calcularmos R̂T R̂ = K − ssT estamos essencialmente
subtraindo de A uma matriz de posto 1, L1 LT1 , o produto externo da primeira
coluna de L pela primeira linha de LT , para então fatorarmos o bloco não nulo
da diferença.
Assim, para completar a fatoração, aplicamos a mesma ideia recursivamente,
ao bloco de dimensão (n − 1) R̂T R̂ = K − ssT e assim por diante, até que o
bloco a ser fatorado seja um escalar positivo.
Cabe destacar que podemos adotar esta abordagem pois
++ ++
A ∈ Sn++ → K ∈ Sn−1 ⇐⇒ R̂T R̂ ∈ Sn−1 .
Exemplo 13 Computar a Fatoração de Cholesky A = LLT = RT R de A dada
abaixo por meio da visão de produtos externos, outer Cholesky.
 
1 −1 3 −4
 −1 5 −1 2 
A=  3 −1 14 −9 


−4 2 −9 22
1. Primeira op, determinamos a primeira linha de LT = R.
√
r11 = a11 = 1
wT
sT = = (−1, 3, −4)
r11

Com isso, calculamos o segundo bloco a ser fatorado RT R:

   
5 −1 2 −1 
R̂T R̂ =  −1

14 −9  −  3  −1 3 −4
2 −9 22 −4
 
4 2 −2
= 2 5 3 
−2 3 6

31
  
4 2 −2
2. Segunda op, desejamos a fatoração de Cholesky de  2 5 3 , ou
−2 3 6
seja, desejamos a segunda linha de LT = R.
√
r22 = 4 = 2
wT
sT = = (1, −1)
2

O termo que sobra RT R para ser fatorado é:

   
T 5 3 1  
R̂ R̂ = − 1 −1
3 6 −1
 
4 4
=
4 5
 
4 4
3. Terceira op, desejamos a fatoração de Cholesky de , ou seja de-
4 5
terminamos a terceira linha de LT = R
√
r33 = 4 = 2
wT
sT = = (2)
2

O termo que sobra RT R:

R̂T R̂ = 5 − (2)(2)T = 1

4. Quarta op, determinamos a quarta linha de LT = R, obtendo a fatoração

de Cholesky da matriz [1].
√
r44 = 1 = 1

Compondo
 as linhaque calculamos, podemos escrever o fator LT = R =
1 −1 3 −4

 2 1 −1  .
 2 2 
1

5 Sistemas lineares malcondicionados

Com as fatorações P A = LU e de Cholesky A = LLT podemos resolver uma
grande variedade de problemas em Álgebra Linear, por exempo, a resolução de
sistemas lineares. Elas são úteis para se resolver boa parte dos sistemas lineares
que usualmente encontramos em aplicações, sobretudo se estruturas de dados
adequadas para representação de matrizes esparsas forem empregadas.

32
 Entretanto, há limitações para o uso destas fatorações. Quando os sistemas
lineares são malcondicionados, outras fatorações devem ser empregadas. Estas
outras fatorações, seja para a resolução de sistemas lineares ou para outros
propósitos, e a base matemática que as fundamenta é o foco do restante do
curso de ALC.
Porém, antes de apresentarmos estas outras fatorações, devemos definir quais
são estes tipos de sistemas lineares que não devemos esperar que as fatorações
básicas sejam capazes de resolver. Devemos caracterizar os sistemas lineares
Ax = b malcondidicionados, que são definidos por matrizes de coeficientes A
malcondicionadas. Na apresentação desta seção, salvo menção contrária, assu-
mimos que A é não singular e que b é diferente de zero. Assim sendo, o sistema
linear Ax = b admite solução única não nula.
O bom ou o mau condicionamento de um sistem linear Ax = b depende
de uma grandeza associada à matriz de coeficientes: o número de condição da
matriz A, κ(A).
Dada uma norma matricial ∥·∥ induzida por uma norma vetorial p, definimos
o número de condição na norma p como

κp (A) = ∥A∥p ∥A−1 ∥p . (12)

Quando o valor de p não for relevante para a análise ou quando for claro
pelo contexto, usaremos κ(A) para fazer menção ao número de condição da ma-
triz. Lembre-se que a norma matricial espectral é induzida pela norma vetorial
p = 2, ou norma Euclideana. Recorde-se também que a norma de Frobenius
é subordinada à norma Euclideana, nas não é induzida por ela. Portanto, não
podemos definir o número de condição de uma matriz, à partir da norma de
Frobenius.
Veja que o número de condição de uma matriz depende da sua inversa,
de forma que o seu cálculo é computacionalmente custoso (no mı́nimo O(n3 )).
Observe também que, pela definição, κ(A) = κ(A−1 ) vale. Além disso, se a
norma matricial espectral for a norma escolhida, a avaliação de κ(A)2 torna-se
mais onerosa. Portanto, é comum utilizamos limites inferiores para κ(A) para
inferir propriedades da matriz. Em particular, sendo A1 , . . . , An as colunas de
A, podemos derivar limites inferiores para o número de condição de A utilizando
qualquer par de ı́ndices de colunas i, j ∈ {1, . . . , n} e computando o lado direito
da seguinte desigualdade:

∥Ai ∥p
κp (A) ≥ , 1 ≤ p ≤ ∞. (13)
∥Aj ∥p
Veja que para i = j, a desigualdade (13) indica que κp (A) ≥ 1, para qualquer
norma matricial induzida por norma vetorial.

Exemplo 14 Duas matrizes notavelmente malcondicionadas são as matrizes

de Hilbert H e de Vandermonde V , indicadas abaixo. O termo geral da entrada
hij da matriz de Hibert é hij = 1/(i + j − 1). Assim, para n = 4, temos
 
1 1/2 1/3 1/4
 1/2 1/3 1/4 1/5 
H4 =   1/3 1/4 1/5 1/6 


1/4 1/5 1/6 1/7

33
 Para definirmos uma matriz de Vandermonde, precisamos definir um vetor
de pontos ou de dados x = (x1 , . . . , xm )T , com entradas distintas par a par:
xi ̸= xj para qualquer par i ̸= j. Já a expressão do termo geral da matriz de V
associada a este conjunto de dados é vij = xj−1
i para i = 1, . . . , m e j = 1, . . . , n.
Veja o caso quadrado:

1 x1 x21 . . . xn−1
 
1
 1 x2 x22 . . . xn−1 
2
V = .
 
.. .. ..
 ..

. . . 
1 xn x2n ... xn−1
n

Note que para qualquer diferença de magnitude entre a mı́nima e a máxima das
entradas em {x1 , . . . , xm }, a potenciação destas entradas aumenta a diferença
de escala. Usando a expressão (13), fica evidente que tanto H quando V são
malcondicionadas.

A utilidade do número de condição, em qualquer norma que seja, consiste

em dizer o quão sensı́vel é a solução x do sistema linear Ax = b, quando as
entradas do sistema linear, A, b ou ambos, são perturbados.
A matriz A é bem condicionada quando seu número de condição, não importa
qual p você use, é pequeno. A definição precisa de pequeno depende de muitos
fatores, por exemplo o número de bits sendo empregado na representação dos
números de ponto flutuante, a quantidade de erros numéricos aceitáveis para
uma dada aplicação e a confiança que temos na qualidade dos dados do problema
(o quão precisos são A, b na modelagem da aplicação representada pelo sistema).
Apesar destas observações, admite-se que um número de condição bom deva
ser inferior a 100. Voltaremos a discutir este aspecto em breve.
O fato é que os dados A, b carregam erros. Seja porque os valores armazena-
dos internamente na máquina são aproximações dos dados verdadeiros (armaze-
namos f l(A), f l(b) e não A, b propriamente), ou porque são dados que vieram do
laboratório, de cálculos computacionais anteriores ou de modelos matemáticos
simplificados para se representar problemas muito complexos. Os dados carre-
gam erros. Assim é razoável pensarmos em A, b como dados verdadeiros de um
sistema linear hipotético, não sujeitos a qualquer tipo dos erros discutidos e em
Ã, b̃ como os dados do sistema linear que efetivamente vamos resolver, pois são
os dados que dispomos para aproximar A, b.
Diante desta perspectiva, informalmente, se b̃ é um vetor muito próximo de
b, as soluções x̃ e x respectivamente de Ãx̃ = b̃ e Ax = b devem ser próximas,
quando o sistema linear é bemcondidiconado. Quando o sistema é malcondici-
onado (isto é, κ(A) é elevado), pequenas variações nos dados do sistema linear
acarretam grandes variações na solução do sistema: ∥x − x̃∥ é muito grande
mesmo quando ∥b − b̃∥ (ou ∥A − Ã∥ ou ambos) é muito pequeno. Esta afirma-
tiva precisa ser demonstrada matematicamente. É o que vamos fazer, mostrar
que κ(A) funciona como uma garantia de que erros muito grandes na resposta
não sejam obtidos se os dados A, b e Ã, b̃ não diferem muito. Para isso, como
a discussão até aqui já sugere, vamos usar normas matriciais e vetoriais para
mensurar a magnitude das perturbações nos dados δ(b) = b − b̃ e δ(A) = A − Ã
e na solução do sistema linear δ(x) = x − x̃.
De inı́cio, vamos considerar o caso em que apenas o vetor b é sujeito à
erros, de forma que A = Ã. Então temos que Ax̃ = b̃ equivale a A(xδ(x)) =

34
b + δ(b). Como Ax = b, o sistema anterior pode ser escrito como Aδ(x) =
δ(b) ou equivalentemente δ(x) = A−1 δ(b). Recordamos que sempre usaremos
normas matriciais induzidas por normais matriciais. Assim sendo, usando as
propriedades de normas matriciais induzidas por normas vetoriais, temos que
∥δ(x)∥ ≤ ∥A−1 ∥|δ(b)∥. Aplicando a mesma relação ao sistema não perturbado
1 1
por erros, b = Ax, temos ∥b∥ ≤ ∥A∥|x∥ → ∥x∥ ≤ ∥A∥ ∥b∥ . Combinando as duas
desigualdades obtemos:

∥δ(x)∥ ∥δ(b)∥
≤ ∥A∥∥A−1 ∥
∥x∥ ∥b∥
∥δ(x)∥ ∥δ(b)∥
≤ κ(A) . (14)
∥x∥ ∥b∥

Veja que a expressão acima mostra que κ(A) é uma grandeza que surge
naturalmente quando tentamos relacionar a perturbação relativa nos dados,
∥δ(b)∥ ∥δ(x)∥
∥b∥ , em função da perturbação observada na solução do sistema linear, ∥x∥ .
Daı́ decorre a definição e utilidade do número de condição.
Observe agora que a desigualdade (14) indica que κ(A) funciona como uma
trava: se κ(A) é pequeno e ∥δ(b)∥ ∥δ(x)∥
∥b∥ também é pequeno, não há como esperar ∥x∥
grande. Por outro lado, mesmo que ∥δ(b)∥ ∥b∥ seja pequeno, se κ(A) é grande, a
perturbação da resposta pode ser grande. E o fato de usarmos normas matriciais
induzidas por normas vetoriais, sempre haverá, para uma matriz A e um b, uma
perturbação δ(b) que faça a desigualdade (14) ser satisfeita de forma justa, na
igualdade (o lado direito e esquerdo da desigualdade assumindo valores iguais).
Usando argumentos semelhantes, podemos demonstrar a validade da se-
guinte desigualdade, caso as perturbações ocorram apenas em A e o vetor b
não seja perturbado:

∥δ(x)∥ ∥δ(A)∥
≤ κ(A)
∥x + δ(x)∥ ∥A∥

Exemplo 15 Neste exemplo, vamos aplicar aplicar uma perturbação pequena

em b e verificar o que acontece com a resposta
 na solução
 do sistema linear 
1000 999 −1 −998 −999
definido por Ax = b onde A = , A = e
999 998 999 −1000
10−3
   
1
b= . Vamos assumir que a perturbação é dada por δb = .
1 0
Para este exemplo, vamos calcular explicitamente o número de condição nas
normas p = 1, ∞. Para tanto, usamos o fato de que ∥A∥1 = ∥A∥∞ = ∥A−1 ∥1 =
∥A−1 ∥∞ = 1999. Logo κ1 (A) = κ∞ (A)  × 10 .
 = 3.992
6

1
A solução de Ax = b é x = . Por outro lado, a solução do
−1
2 × 10−3
 
sistema linear perturbado Ax = (b + δb ) é x̂ = . Veja que
  −1 × 10−3
−0.998
δx = . Portanto ∥δb∥1 = 1.0 × 10−3 , ∥δx∥1 = 1.997, ∥δx∥ 1
∥δb∥1 =
0.999
1.997 × 10+3 . A perturbação na resposta é cerca de 2000 vezes maior que a
perturbação nos dados.

35
 A matriz A do exemplo não é singular, pois det(A) = 1. Porém, verificamos
também que as duas restrições do sistema linear são praticamente linearmente
dependentes. E esse fato nos motiva a produzir uma discussão adicional.
Observe que o número de condição κ(A) é uma propriedade da matriz e
nada tem a ver com a precisão da máquina onde eventualmente resolveremos
a solução do sistema linear. É uma propriedade inexorável da matriz. É dela,
não depende da máquina.
Então por que afirmamos que sistemas lineares malcondicionados são difı́ceis
de serem resolvidos ou de serem fatorados ? Se usarmos aritmética de precisão
infinita, um sistema linear malcondicionado (em que A é inversı́vel) não é pior
do que outro bem condicionado. Porém, no processo de fatoração das ma-
trizes, inevitavelmente calculamos fatores sujeitos a erros numéricos. E estas
perturbações levam a perturbações grandes nas respostas.
No caso de matrizes malcondicionadas, esses erros numéricos são grandes
o suficiente para, muitas vezes, introduzir dependência linear onde não há,
matematicamente, de forma exata, depedência linear. Para ser mais preciso,
recorde-se do exemplo (9) que apresentamos. A matriz A daquele exemplo é
bem condicionada, κ(A)1 = 4. Naquele exemplo, empregamos a fatoração LU
sem troca de linhas e diante disso, o fator Ũ que foi obtido é uma matriz onde as
duas linhas são praticamente linearmente dependentes, quando a matriz A origi-
nal é uma matriz em que suas linhas são claramente linearmente independentes.
E em função disso, sugerimos a troca de linhas, que para matrizes bem condici-
onadas é uma ideia numérica eficaz para produzir soluções com boa qualidade
numérica.
Porém, o potencial de redução de erros numéricos da permutação de linhas
quando a matriz é malcondicionada é limitado. Para uma matriz malcondi-
cionada, mesmo a introdução de mecanismos de pivoteamento resulta na in-
trodução de dependência linear no sistema equivalente U x = y a ser resolvido.
Isso porque as entradas das matrizes Aj−1 : j = 1, . . . , n − 1 que competem
para definir a linha pivotal já são substancialmente diferentes daquelas que se-
riam obtidas se aritmética de precisão infinita tivesse sido empregada. Isso faz
com que mecanismos distintos da Eliminação de Gauss sejam empregados. Os
mecanismos empregados na Eliminação de Gauss para produzir um sistema tri-
angular superior U x = y equivalente a Ax = b, qual seja, empregar combinações
lineares das linhas do sistema e troca de linhas, são propensos a erros.

5.1 Resolvendo um sistema linear definido por uma matriz

de Vandermonde usando P A = LU .
O nosso objetivo nesta seção é ilustrar o efeito do malcondicionamento da ma-
triz na qualidade da solução numérica produzida pela fatoração P A = LU para
se resolver o sistema linear. Para tanto, construiremos um sistema linear mal-
condicionado, definido por uma matriz de Vandermonde. O sistema linear será
construı́do partindo da solução desejada para o sistema linear, um vetor n-
dimensional de 1’s. Resolveremos o sistema linear para diversos valores de n e
verificaremos a diferença entre a resposta numérica e a resposta que esperávamos
obter.
Para criar este experimento, vamos considerar o polinômio de grau n − 1 na

36
variável t:
n−1
X
p(t) = ai ti
i=0

onde a0 = a1 = · · · = an−1 são os coeficientes do polinômio. No nosso experi-

mento, vamos fixar todos os coeficientes em 1. Isto é, o polinômio que vamos
considerar é:

p(t) = 1 + t + t2 + · · · + tn−1 (15)

Agora, vamos arbitrariamente escolher um conjunto de n (o número de coe-
ficientes do polinômio de grau n−1) abiscissas t e, à partir delas, vamos calcular
suas correspondentes ordenadas p(t). Com as n tuplas (t, p(t)) construiremos um
sistema linear no qual a matriz de coeficientes é uma matriz de Vandermonde, o
vetor de termos independentes consiste no vetor formando pelos valores assumi-
dos pelo polinômio para cada uma destas entradas e o vetor x que procuramos é
o vetor dos coeficientes de um polinômio que interpola todos estes pontos. Veja
que já sabemos a resposta correta para o sistema linear: o vetor n-dimensional
(1, . . . , 1)T .
A criação destes dados seguirá os seguintes passos:

• Parametrizamos t = i + 1, para diversos valores distintos de i.

• Para cada valor de t = i + 1, p(t) pode ser reescrito como o valor da soma
de uma Progressão Geométrica, cuja expressão analı́tica é desenvolvida
abaixo:
n−1
X
p(i + 1) = (i + 1)0 + (i + 1) + (i + 1)2 + · · · + (i + 1)n−1 = (i + 1)j
j=0

Ou seja, para um dado valor de i, p(i + 1) corresponde à soma dos termos

de uma Progressão Geométrica de n termos, com o primeiro termo igual
a 1 e razão (i + 1). Então podemos escrever:

n−1
X (1 + i)n − 1
p(i + 1) = (i + 1)j =
j=0
i

A Figura 7 apresenta o código da geração do sistema linear que desejamos

resolver (função GeraVandermonde) e o procedimento GeraExperimento, que
seá executado com o seguinte vetor de dados de entrada para valores de n:
(5, 8, 10, 15, 20)T . Observe que para cada n ∈ valoresden, o procedimento cal-
cula ∥x∥∞ , a norma do vetor x, solução numérica do sistema linear. Veja na
Figura 8, os valores das normas das soluções x encontradas. Para n = 10 em
diante, em nada estas normas conferem com a norma infinito de um vetor de
uns, que é 1. Em particular, para valores de n = 15, 20, os fatores encontrados
na fatoração LU (usando a implementação profissional disponı́vel no scilab)
são absolutamente distintos dos valores corretos. Isso ocorre pois a matriz de
coeficientes é extremamente mal condicionada.
O problema da fatoração LU para lidar com matrizes malcondicionadas está
na ideia central do método, que é construir combinaçõesl lineares das linhas de

37
function [A,b] = GeraVandermonde(n)
for i = 1:n
for j = 0:n-1
A(i,j+1) = (1 + i)^j;
end
b(i) = ((1 + i)^n - 1)/i;
end
endfunction

function [normas] = GeraExperimento(valoresn)

s = size(valoresn)
for i = 1:s(1)
[A,b] = GeraVandermonde(valoresn(i))
[L,U,P] = lu(A)
[m,n] = size(A)
[y] = SubsSucessivas(L,P*b,n)
[x] = SubsRetroativas(U,y,n)
printf("n = %d %8.7E \n",valoresn(i),norm(x,’inf’))
normas(i) = norm(x,%inf)
end
endfunction

Figura 7: Procedimentos para o experimento numérico com a matriz de Van-

dermonde.

Ax = b, visando triangularizar a matriz. Recorde-se das transformações lineares

Mn−1 Mn−2 . . . M1 A = U . Estas são instáveis, pois as matrizes de multiplicado-
res Mj possuem entradas obtidas por divisões pelo elemento pivot, que pode ter
magnitude muito pequena, ainda que usemos o pivoteamento de colunas. Assim,
em boa parte do restante deste curso, apresentaremos outras ideias para pro-
duzirmos outras fatoração matriciais, mais estáveis numericamente, no espı́rito
que narramos aqui.

38
-->valoresden’
ans =
5. 8. 10. 15. 20.
-->normas = GeraExperimento(valoresden)
n = 5 1.0000000E+00
n = 8 1.0000001E+00
n = 10 1.0003594E+00
n = 15 5.3951429E+05
n = 20 2.0111339E+18
normas =
1.
1.0000001
1.0003594
539514.29
2.011D+18

Figura 8: Solução obtida para o experimento da Figura 7.

39