Vectores y Matrices
Curso 2017-18
Notación
a11 a12 . . . a1n
a21 a22 . . . a2n
A= . .. .. .. , A = [aij ] 1≤i≤m
.. . . . 1≤j≤n
am1 am2 . . . amn
aij = elemento (i, j) de A.
Tamaño u orden de A = m × n. Si m = n, A cuadrada.
i-sima fila y j-sima columna:
a1j
a2j
i
a = [ai1 ai2 . . . ain ] y aj = . ,
..
amj
Notación de MATLAB: A(i,:) y A(:,j).
En general aij ∈ R, anillo conmutativo con elemento 1.
Rm×n = cto de las matrices de tamaño m × n con elementos en R.
2
Operaciones con matrices
Si A = [aij ] ∈ Rm×n y B = [bij ] ∈ Rm×n :
A + B = [aij + bij ] ∈ Rm×n .MATLAB:A+B
y si A = [aij ] ∈ Rm×n y B = [bij ] ∈ Rn×p , entonces
" n #
X
AB = aik bkj ∈ Rm×p .MATLAB:A*B
k=1
(Rn×n , +, ·) es un anillo no conmutativo con identidad
1 0 ... 0
0 1 ... 0
matriz identidad de orden n.
In = . .. . . .. = [δij ],
.. . . . MATLAB:eye(n)
0 0 ... 1
Si x ∈ R y A = [aij ] ∈ Rm×n , xA = [xaij ] ∈ Rm×n .
MATLAB: x.∗A o x∗A
3
Algunos tipos de matrices
Si A = [aij ]
AT = [aji ] transpuesta de A. MATLAB: A’ si R = R y A.’ si R = C.
Si R = C, A = [aij ] conjugada de A. (z = x + iy ⇒ z = x − iy )
T
Si R = C, A∗ = A . MATLAB: A’
1 2 3
i 2i
Simétricas: A = AT : 2 0 4 ,
2i 4
3 4 −1
2 −i 2 + i
Hermı́ticas: A = A∗ : i 3 5
2−i 5 1
0 1 + 2i
Skew-simétricas: AT = −A:
−1 − 2i 0
Skew-hermı́ticas: A∗ = −A.
4
Más tipos de matrices
d1 0 · · · 0
0 d2 · · · 0
Diagonales: . .. . . .. = Diag(d1 , . . . , dn )
.. . . .
0 0 ··· dn
Matlab: diag
A1 0 . . . 0
0 A2 . . . 0
Diagonal por bloques: .. .. . . . .
. . . ..
0 0 . . . Ap
Matlab: blkdiag
a11 a12 · · · a1n
0 a22 · · · a2n
Triangular superior: . .
..
.
.. . .. ..
.
0 0 · · · ann
Matlab: triu
Triangular inferior, por bloques, etc.
5
Matlab: tril
Matrices de permutación
P ∈ Rn×n matriz de permutación si en cada fila y columna hay un
elemento igual a 1 y todos los demás son cero: Permutan filas y columnas.
0 0 1 0 1 0
P = 1 0 0 y PT = 0 0 1
0 1 0 1 0 0
1 3
P 2 = 1 y 1 2 3 PT = 3 1 2 .
3 2
Si σ = (i1 , i2 , . . . , in ) ∈ Sn ; i. e. σ(k) = ik , Pσ = [δiσ(j) ]= matriz con un 1
en (σ(i), i) y todo lo demás cero . En el ejemplo, σ = (2, 3, 1).
P T = Pσ−1 = P −1
Pσ A = [aσ−1 (i)j ] y AP T = [aiσ−1 (j) ].
Matlab: P=eye(n); P=P(:,σ)
6
Submatrices
Si A = [aij ] ∈ Rm×n , B = [bij ] ∈ Rp×q submatriz de A si existen ı́ndices
(i1 , . . . , ip ) y (j1 , . . . , jq ) tales que 1 ≤ i1 < · · · < ip ≤ m, 1 ≤ j1 < · · · < jq ≤ n y
air js = brs , 1 ≤ r ≤ p, 1 ≤ s ≤ q.
rb11 rb12 rb13 i1
rb21 rb22 rb23 i2
A=
rb31 rb32 rb33 i3
rb41 rb42 rb43 i4
j1 j2 j3
Usaremos notación de MATLAB: si u = (i1 , . . . , ip ) y v = (j1 , . . . , jq ):
B=A(u,v)
Hay configuraciones que no son submatrices. 7
Matrices y aplicaciones lineales
R=F
f : V1 → V2
Fijadas bases B1 = {v1 , . . . , vn }, B2 = {u1 , . . . , um }
m
X
f (vj ) = aij ui , j = 1, . . . , n
i=1
determina A = [aij ] ∈ Fm×n : matriz de f respecto de B1 y B2 . Y
recı́procamente.
A ∈ Fm×n puede verse como una aplicación lineal entre los espacios
vectoriales Fn y Fm :
A : Fn → Fm
x ; Ax
La matriz de esta aplicación lineal es A cuando en Fm y Fn se toman las
bases canónicas.
8
Rango y nulidad
A ∈ Fm×n
Im(A) = {y ∈ Fm |y = Ax para algún x ∈ Fn }
Ker(A) = {x ∈ Fn |Ax = 0}
rang(A) = dim Im A MATLAB: rank
ν(A)= nulidad de A= dim Ker A
Primer teorema de isomorfı́a
Fn
Im A ∼
= ⇒ ν(A) = n − rang(A).
Ker A
Desigualdad de Sylvester (1884)
A ∈ Fm×n , B ∈ Fn×p
rang(A) + rang(B) − n ≤ rang(AB) ≤ mı́n(rang(A), rang(B))
Matrices de rango completo
A tiene rango completo ⇔ rang(A) = mı́n{m, n}.
¿Qué significa que A es de rango completo?
m ≥ n ⇒ dim Im A = n ⇔ Ker A = {0} ⇔ A es inyectiva (invertible
por la izquierda: BA = In ).
m ≤ n ⇒ dim Im A = m ⇔ Im A = Fm ⇔ A es suprayectiva
(invertible por la derecha: AB = Im ).
m = n ⇒ A biyectiva (invertible: AB = BA = In ).
10
Matrices invertibles
Para A ∈ Fn×n (F cuerpo) las siguientes condiciones son equivalentes:
1 A tiene inversa A−1 . (MATLAB: inv)
2 rang(A) = n.
3 Im A = Fn .
4 ν(A) = 0.
5 Ker A = {0}.
6 0 no es un valor propio (autovalor) de A. (MATLAB: eig)
7 det(A) 6= 0. (MATLAB: det)
A no singular: det(A) 6= 0.
A no singular ⇔ A invertible (F cuerpo)
λ ∈ C valor propio de A ∈ Fn×n :
∃x ∈ Cn (x 6= 0) tal que Ax = λx
x= vector propio de A asociado a λ.
11
¿Qué significa Ax = b?
b es la imagen de x por A
b es el resultado de multiplicar A por x.
Xn
bi = aij xj , i = 1, . . . , m.
j=1
b1 a11 a12 a1n
b2 a21 a22 a2n
.. = .. x1 + .. x2 + · · · + .. xn .
. . . .
bm am1 am2 amn
b = x1 a1 + x2 a2 + · · · + xn an .
b es una combinación lineal de las columnas de A cuyos
coeficientes son las componentes del vector x.
n
X
b ∈ Im A ⇔ b = xi ai para algún vector x
i=1
12
La Imagen y el Span
Notación: Si a1 , . . . , an ∈ Fm
( n )
X
< a1 , . . . , an >= Span (a1 , . . . , an ) = xi ai |xi ∈ F .
i=1
de F generado por las columnas de A; i.e., si
m
es el subespacio
Im(A)
A = a1 a2 · · · an , Im A =< a1 , . . . , an >.
rang(A) = dim < a1 , . . . , an >
Ker A está formado por los vectores cuyas componentes son los
coeficientes del vector 0 como combinación lineal de las columnas de A:
0 = x1 a1 + x2 a2 + · · · + xn an .
13
Producto de matrices
¿ Qué significa C = AB?
B A
C es la matriz de la composición Rp −→ Rn −→ Rm
C es el resultado de multiplicar A por B.
cij = ai1 b1j + ai2 b2j + · · · + ain bnj ,
c1j a11 a12 ··· a1n b1j
c2j a21 a22 ··· a2n b2j
cj = . = . .. .. .. .. = Abj , j = 1, . . . , p.
.. .. . . . .
cmj am1 am2 · · · amn bnj
c1 c2 ··· cp = Ab1 Ab2 ··· Abp ,
c i = ai B.
A actúa sobre B: combinaciones lineales en las filas de B, y
B actúa sobre A: combinaciones lineales en las columnas de A. 14
Producto interno y externo de vectores
Si u, v ∈ Rn :
Producto interno: u T v
v1
T
v2 X
n
u v = u1 u2 · · · un . = ui vi ∈ R.
..
i=1
vn
producto externo: uv T ∈ Rn×n
u1 u1 v T
u2 T
T u2 v
uv = . v1 v2 · · · vn = . = v1 u v2 u · · · vn u .
.. ..
un un v T
uv T es una matriz de rango 1 y todas las matrices
de rango 1 son ası́.
15
Otra vez el producto de matrices
P
n
(AB)ij = aik bkj = ai bj = (ai0 )T bj , ai0 es la i-ésima columna de AT
k=1
y bj la j-ésima columna de B.
(ak b k )ij = aik bkj
AB = a1 (b10 )T + a2 (b20 )T + · · · + an (bn0 )T ,
16
Matrices elementales. Tipo I
i j
↓ ↓
1
..
.
1 ... α
← i
.. ..
E1 (α) = In + Eij (α) = . .
1 ← j
..
.
1
a1
..
.
i
a + αaj
..
E1 (α)A = . , AE1 (α) = a1 ··· ai ··· aj + αai ··· an .
aj
..
.
am 17
Matrices elementales. Tipo II
i
↓
1
..
.
1
E2 (α) =
α
← i
1
..
.
1
a1
..
.
i
E2 (α)A =
αa
y AE2 (α) = a1 ··· αai ··· an .
..
.
am
18
Matrices elementales. Transposiciones
i j
↓ ↓
1
..
.
1
0 ... 1
← i
.. .. ..
P = . . .
1 ... 0 ← j
1
..
.
1
1
a
..
.
j
a
PA = ... y AP = a1 ··· aj ··· ai ··· an .
ai
.
..
am 19
Eliminación Gaussiana (EG)
Ax = b, A ∈ Fn×n , b ∈ Fn×1
A b = A(1) b (1) =⇒ A(n) b (n)
A(n) = U triangular superior.
n − 1 etapas y en cada etapa k = 1 : n − 1
" #
(k) (k)
A11 A12 (k)
A(k) = (k) , A11 ∈ F
(k−1)×(k−1)
triangular superior
0 A22
Objetivo en la etapa k: sustituir por ceros los elementos de
A(k) (k + 1 : n, k) restando a la fila i = k + 1 : n la fila k multiplicada por
(k)
aik
mik = (k) (MULTIPLICADOR):
akk
(k+1) (k) (k)
aij = aij − mik akj , i, j = k + 1, . . . , n.
(k+1) (k) (k)
bi = bi − mik bk , i = k + 1, . . . , n.
20
Factorización
LU
T
Si mk = 0 · · · 0 mk+1k ··· mnk y Mk = In − mk ekT :
Mk A(k) = A(k+1)
Entonces
Mn−1 Mn−2 . . . M2 M1 A = A(n) = U triangular superior
A = M1−1 M2−1 . . . Mn−1
−1
U, Mk−1 = In + mk ekT
1
m21 1
.
. . .
−1 −1 −1
L := M1 M2 . . . Mn−1 = . m32 .
.. .. . ..
. .
mn1 mn2 · · · mnn−1 1
A = LU Factorización LU de A
21
Algortimo LU
Factorización LU
Dato: A ∈ Fn×n
Objetivo: calcular factorizarión LU de A.
• L =eye(n)
• for k = 1 : n − 1
• for i = k + 1 : n
• lik = aik /akk
• aik = 0
• for j = k + 1 : n
• aij = aij − lik akj
• end for
• end for
• end for
•U=A
Coste: ∼ 23 n3
22
Pivoteo parcial por columnas
al comienzo de la etapa k, se intercambian las filas k y r estando r
determinada por la condición:
(k) (k)
|ark | = máx |aik |
k≤i≤n
Mn−1 Pn−1 Mn−2 Pn−2 . . . M2 P2 M1 P1 A = U,
U = M4 P4 M3 P3 M2 P2 M1 P1 A
= M4 · P4 M3 P4 · P4 P3 M2 P3 P4 · P4 P3 P2 M1 P2 P3 P4 · P4 P3 P2 P1 ·A
| {z } | {z } | {z } | {z }
M30 M20 M10 P
= M4 M30 M20 M10 PA
Mk0 = Pn−1 Pn−2 . . . Pk+1 (In − mk ekT )Pk+1 . . . Pn−2 Pn−1
= In − Pn−1 Pn−2 . . . Pk+1 mk ekT
= In − mk0 ekT ,
1
m21 0
1
.
.. ..
PA = LU, L = M10−1 M20−1 . . . Mn−2
0−1 −1
Mn−1 = m 0
32
.
.. .. ..
. . .
0 0 23
mn1 mn2 · · · mnn−1 1
Algoritmo Factorización LU con pivoteo. Versión 1
Factorización LU con pivoteo parcial. Versión 1
Dato: A ∈ Fn×n , det A 6= 0
Objetivo: calcular P,L U tales que PA = LU.
• L =eye(n)
• for k = 1 : n − 1
• Calcúlese el máx |A(k : n, k)| y su posición, r .
• Permútense las filas k y r de A, L y P
for i = k + 1 : n
lik = aik /akk
aik = 0
for j = k + 1 : n
aij = aij − lik akj
end for
end for
end for
U=A 24
Algoritmo Factorización LU (LUTXFOR). Versión 1
function [L,U,P] = lutxforv1(A)
%LUTXFORV1, para una matriz invertible dada A, devuelve matrices L,
%triangular inferior con 1 en la diagonal,
%U, triangular superior y P de permutación
%tal que PA=LU. En vez trabajar con la matriz P se trabaja
%con un vector p que almacena las transposiciones y a partir
%del cual se construye P. La matriz L se define
%y actualiza de forma explı́cita y se utilizan bucles FOR para
%todos los contadores.
% Calculamos el tama~
no de A
[n,n] = size(A);
%Iniciamos p=[1 2 ... n]^T
p = (1:n)’;
% y L
L=eye(n);
%Transformaciones elementales en las columnas 1, 2,..., n-1.
for k = 1:n-1
%Calculamos el máximo, m, (en valor absoluto) de los elementos
% A(k,k), A(k+1,k),...,A(n,k), y la posición, r, donde se encuentra
[m,r] = max(abs(A(k:n,k)));
% Calculamos la posición del máximo en toda la columna k-ésima de A
r = r+k-1;
%Si el máximo fuera cero es porque toda la columna es cero y A no serı́a
%invertible
if (A(r,k) ~= 0)
%Si el máximo no está en la diagonal permutamos
if (r ~= k)
%permutamos las filas k y r en A, L y p
A([k r],k:n) = A([r k],k:n);
L([k r],1:k-1) = L([r k],1:k-1);
p([k r]) = p([r k]); 25
end
Algoritmo Factorización LU (LUTXFOR). Versión 1. Cont.
%L está fromada por los multiplicadores
for i = k+1:n;
L(i,k) = A(i,k)/A(k,k);
% Por debajo de (k,k) los elementos de A en la columna k serán 0
A(i,k)=0;
% Realizamos sobre las columnas j=k+1,...,n la transformación elemental
for j = k+1:n;
A(i,j) = A(i,j) - L(i,k)*A(k,j);
end
end
end
end
%A se ha convertido en triangular superior: es la U
U=A;
% Construı́mos P a partir de p: las filas 1:n de P son las de la permutación
P=eye(n); P=P(p,:);
26
Algoritmo Factorización LU (LUTXFOR). Final
function [L,U,P] = lutxfor(A)
[n,n] = size(A);
p = (1:n)’;
for k = 1:n-1
[m,r] = max(abs(A(k:n,k)));
r = r+k-1;
if (A(r,k) ~= 0)
if (r ~= k)
A([k r],:) = A([r k],:);
p([k r]) = p([r k]);
end
for i = k+1:n;
A(i,k)=A(i,k)/A(k,k);
for j = k+1:n;
A(i,j) = A(i,j) - A(i,k)*A(k,j);
end
end
end
end
%En la parte triangular inferior de A (sin la diagonal) está L
L = tril(A,-1) + eye(n,n);
%Y en la parte triangular superior (incluyendo la diagonal) está U
U = triu(A);
P=eye(n); P=P(p,:);
27
El algoritmo de la factorización LU
function [L,U,p] = lutx(A)
[n,n] = size(A);
p = 1:n;
for k = 1:n-1
[m,r] = max(abs(A(k:n,k)));
r = r+k-1;
if (A(r,k) ∼= 0)
if (r ∼= k)
A([k r],:) = A([r k],:);
p([k r]) = p([r k]);
end
i = k+1:n;
A(i,k) = A(i,k)/A(k,k);
j = k+1:n;
A(i,j) = A(i,j) - A(i,k)*A(k,j);
end
end
L = tril(A,-1) + eye(n,n);
U = triu(A);
28
Jerarquı́a en la memoria de los ordenadores y BLAS
Registros
Caché
RAM
Disco
Basic Linear Algebra Subroutines (BLAS)
Operación Definición flops (f ) memory ref. (m) Ratio q = f /m
BLAS1 y = αx + y 2n 3n + 1 2/3
BLAS2 y = Ax + y 2n2 n2 + 3n 2
BLAS3 C = AB + C 2n3 4n2 n/2
BLAS1 (saxpy): productos internos, escalar por vector, . . .
BLAS2: matriz por vector, sistemas triangulares, A + xy t ,. . .
BLAS3: matriz por matriz, sistemas triangulares con muchos vectores
independientes, . . .
m tmem 1 tmem 29
f · tarit + m · tmem = f · tarit · 1 + f tarit = f · tarit · 1 + q tarit
