Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.

Álgebra lineal

Matrices en el entorno MATLAB

Este tema contiene una introducción a la creación de matrices y la realización de cálculos de matrices básicos en MATLAB®.

El entorno MATLAB utiliza el término matriz para indicar una variable que contiene números reales o complejos dispuestos en una cuadrícula bidimensional. Un arreglo es, más generalmente, un vector, una matriz o una cuadrícula de números de más dimensiones. Todos los arreglos de MATLAB son rectangulares. Esto quiere decir que los vectores que lo componen a lo largo de cualquier dimensión tienen todos la misma longitud. Las operaciones matemáticas que se definen en matrices conforman la temática del álgebra lineal.

Creación de matrices

MATLAB tiene múltiples funciones que crean diferentes tipos de matrices. Por ejemplo, puede crear una matriz simétrica con entradas basadas en el triángulo de Pascal:

A = pascal(3)
A =
       1     1     1
       1     2     3
       1     3     6

O bien, puede crear una matriz del cuadrado mágico asimétrica, cuyas filas y columnas sumen lo mismo:

B = magic(3)
B =
       8     1     6
       3     5     7
       4     9     2

Otro ejemplo es una matriz rectangular de 3 por 2 de números enteros aleatorios. En este caso, la primera entrada en randi describe el rango de posibles valores para los enteros, y las segundas dos entradas describen la cantidad de filas y columnas.

C = randi(10,3,2)
C =

     9    10
    10     7
     2     1

Un vector columna es una matriz de m por 1, un vector fila es una matriz de 1 por n y un escalar es una matriz de 1 por 1. Para definir una matriz de forma manual, use corchetes [ ] para indicar el inicio y el final del arreglo. Dentro de los corchetes, use un punto y coma ; para indicar el final de una fila. En el caso de un escalar (matriz de 1 por 1), los corchetes no son necesarios. Por ejemplo, estas instrucciones producen un vector columna, un vector fila y un escalar:

u = [3; 1; 4]

v = [2 0 -1]

s = 7
u =
       3
       1
       4

v =
       2     0    -1

s =
       7

Para obtener más información sobre la creación y el trabajo con matrices, consulte Crear, concatenar y expandir matrices.

Adición y sustracción de matrices

La adición o la sustracción de matrices y arreglos se lleva a cabo elemento por elemento. Por ejemplo, si se suma A y B, y luego se sustrae A del resultado, se recupera B:

X = A + B
X =
       9     2     7
       4     7    10
       5    12     8
Y = X - A
Y =
       8     1     6
       3     5     7
       4     9     2

Tanto la adición como la sustracción requieren que ambas matrices tengan dimensiones compatibles. Si las dimensiones son incompatibles, se produce un error:

X = A + C
Error using  + 
Matrix dimensions must agree.

Para obtener más información, consulte Matriz frente a operaciones matriciales.

Productos y trasposición de vectores

Un vector fila y un vector columna de la misma longitud se pueden multiplicar en cualquier orden. El resultado es un escalar, conocido como producto interno, o una matriz, conocida como producto externo:

u = [3; 1; 4];
v = [2 0 -1];
x = v*u
x =

     2
X = u*v
X =

     6     0    -3
     2     0    -1
     8     0    -4

En el caso de matrices reales, la operación de trasposición intercambia aij y aji. En el caso de matrices complejas, también se debe decidir si tomar o no el conjugado complejo de entradas complejas del arreglo para formar la trasposición conjugada compleja. MATLAB utiliza el operador apóstrofo (') para realizar una trasposición conjugada compleja, y el operador punto y apóstrofo (.') para trasponer sin conjugación. En el caso de las matrices que contienen solo elementos reales, los dos operadores devuelven el mismo resultado.

La matriz de ejemplo A = pascal(3) es simétrica y por lo tanto A' es igual a A. Sin embargo, B = magic(3) no es simétrica, por lo que los elementos de B' se reflejan a lo largo de la diagonal principal:

B = magic(3)
B =

     8     1     6
     3     5     7
     4     9     2
X = B'
X =

     8     3     4
     1     5     9
     6     7     2

En el caso de los vectores, la trasposición convierte un vector fila en un vector columna (y viceversa):

x = v'

x =
       2
       0
      -1

Si tanto x como y son vectores columna reales, entonces, el producto x*y no está definido, pero los dos productos

x'*y

y

y'*x

producen el mismo resultado escalar. Esta cantidad se utiliza con tanta frecuencia que tiene tres nombres diferentes: producto interno, producto escalar o producto punto. Incluso existe una función dedicada para productos punto llamada dot.

Para obtener un vector o matriz complejo, z, la cantidad z' no solo traspone el vector o la matriz, sino que también convierte cada elemento complejo en su conjugado complejo. Es decir, el signo de la parte imaginaria de cada elemento complejo cambia. Por ejemplo, considere la matriz compleja

z = [1+2i 7-3i 3+4i; 6-2i 9i 4+7i]
z =

   1.0000 + 2.0000i   7.0000 - 3.0000i   3.0000 + 4.0000i
   6.0000 - 2.0000i   0.0000 + 9.0000i   4.0000 + 7.0000i

La trasposición conjugada compleja de z es:

z'
ans =

   1.0000 - 2.0000i   6.0000 + 2.0000i
   7.0000 + 3.0000i   0.0000 - 9.0000i
   3.0000 - 4.0000i   4.0000 - 7.0000i

La trasposición compleja no conjugada, donde la parte imaginaria de cada elemento conserva su signo, se indica mediante z.':

z.'
ans =

   1.0000 + 2.0000i   6.0000 - 2.0000i
   7.0000 - 3.0000i   0.0000 + 9.0000i
   3.0000 + 4.0000i   4.0000 + 7.0000i

Para los vectores complejos, el producto escalar x'*y y el producto escalar y'*x son conjugados complejos entre sí. El producto escalar x'*x de un vector complejo consigo mismo es real.

Multiplicación de matrices

La multiplicación de matrices se define de una forma que refleja la composición de las transformaciones lineales subyacentes y permite una representación compacta de sistemas de ecuaciones lineales simultáneas. El producto matricial C = AB está definido cuando la dimensión de columna de A es igual a la dimensión de fila de B, o cuando una de ellas es un escalar. Si A es m por p y B es p por n, su producto C es m por n. En realidad, el producto puede definirse a través de bucles for, notación de colon (dos puntos) y productos punto de vectores de MATLAB:

A = pascal(3);
B = magic(3);
m = 3; 
n = 3;
for i = 1:m
     for j = 1:n
        C(i,j) = A(i,:)*B(:,j);
     end
end

MATLAB utiliza un asterisco para expresar la multiplicación de matrices, como en C = A*B. La multiplicación de matrices no es conmutativa, es decir, comúnmente, A*B no es igual a B*A:

X = A*B
X =
      15    15    15
      26    38    26
      41    70    39
Y = B*A
Y =
      15    28    47
      15    34    60
      15    28    43

Una matriz se puede multiplicar a la derecha por un vector columna y a la izquierda por un vector fila:

u = [3; 1; 4];
x = A*u
x =

     8
    17
    30
v = [2 0 -1];
y = v*B
y =

    12    -7    10

Las multiplicaciones de matrices rectangulares deben satisfacer las condiciones de compatibilidad dimensional: En vista de que A es de 3 por 3 y C es de 3 por 2, puede multiplicarlas para obtener un resultado de 3 por 2 (la dimensión interior común se cancela):

X = A*C
X =

    24    17
    47    42
    79    77

Sin embargo, la multiplicación no funciona en el orden inverso:

Y = C*A
Error using  * 
Incorrect dimensions for matrix multiplication. Check that the number of columns 
in the first matrix matches the number of rows in the second matrix. To perform 
elementwise multiplication, use '.*'.

Puede multiplicar cualquier cosa con un escalar:

s = 10;
w = s*y
w =

   120   -70   100

Cuando multiplica un arreglo por un escalar, el escalar implícitamente se expande para adoptar el mismo tamaño que la otra entrada. A menudo, esto se conoce como expansión escalar.

Matriz identidad

La notación matemática generalmente aceptada utiliza la letra mayúscula I para nombrar a las matrices identidad, es decir, matrices de diversos tamaños con unos en la diagonal principal y ceros en los demás lugares. Estas matrices tienen la propiedad de que AI = A y IA = A siempre que las dimensiones sean compatibles.

La versión original de MATLAB no podía usar I para este fin, ya que no distinguía entre mayúsculas y minúsculas. Además, i ya servía como subíndice y como la unidad imaginaria. Por lo tanto, se introdujo un juego de palabras del idioma inglés. La función

eye(m,n)

devuelve una matriz identidad rectangular de m por n, mientras que eye(n) devuelve una matriz identidad cuadrada de n por n.

Inversa de la matriz

Si una matriz A es cuadrada y no singular (determinante distinto de cero), entonces, las ecuaciones AX = I y XA = I tienen la misma solución X. Esta solución se denomina la inversa de A y se indica mediante A-1. Tanto la función inv como la expresión A^-1 calculan la inversa de la matriz.

A = pascal(3)
A =
       1     1     1
       1     2     3
       1     3     6
X = inv(A)
X =

    3.0000   -3.0000    1.0000
   -3.0000    5.0000   -2.0000
    1.0000   -2.0000    1.0000
A*X
ans =

    1.0000         0         0
    0.0000    1.0000   -0.0000
   -0.0000    0.0000    1.0000

El determinante calculado por la función det es una medida del factor de escalado de la transformación lineal descrita por la matriz. Cuando el determinante es exactamente cero, la matriz es singular y no existe la inversa.

d = det(A)
d =

     1

Algunas matrices son casi singulares y, pese al hecho de que existe una matriz inversa, el cálculo es susceptible de devolver errores numéricos. La función cond calcula el número de condición para la inversión, lo que indica el grado de precisión de los resultados de la inversión de la matriz. El número de condición oscila desde 1 para una matriz numéricamente estable hasta Inf para una matriz singular.

c = cond(A)
c =

   61.9839

Rara vez es necesario formar la inversa explícita de una matriz. A la hora de resolver el sistema de ecuaciones lineales Ax = b, con frecuencia se utiliza erróneamente la función inv. La mejor forma de resolver esta ecuación, desde el punto de vista del tiempo de ejecución y de la precisión numérica, es usando el operador de barra invertida x = A\b de la matriz. Para obtener más información, consulte mldivide.

Producto tensorial de Kronecker

El producto de Kronecker, kron(X,Y), de dos matrices es la matriz de mayor tamaño que se obtiene a partir de todos los productos posibles de los elementos de X con los elementos de Y. Si la matriz X es m por n, y la matriz Y es p por q, entonces kron(X,Y) es mp por nq. Los elementos están dispuestos de manera tal que cada elemento de X se multiplica por la matriz Y completa:

[X(1,1)*Y  X(1,2)*Y  . . .  X(1,n)*Y
                     . . .
 X(m,1)*Y  X(m,2)*Y  . . .  X(m,n)*Y]

El producto de Kronecker a menudo se utiliza con matrices de ceros y unos para crear copias repetidas de matrices pequeñas. Por ejemplo, si X es la matriz de 2 por 2

X = [1   2
     3   4]

e I = eye(2,2) es la matriz identidad de 2 por 2, entonces:

kron(X,I)
ans =

     1     0     2     0
     0     1     0     2
     3     0     4     0
     0     3     0     4

y

kron(I,X)
ans =

     1     2     0     0
     3     4     0     0
     0     0     1     2
     0     0     3     4

Además de la función kron, también sirven para replicar arreglos las funciones repmat, repelem y blkdiag.

Normas de vectores y matrices

La norma-p de un vector x,

xp=(|xi|p)1p,

se calcula mediante norm(x,p). La operación está definida para cualquier valor de p > 1, pero los valores más comunes de p son 1, 2 e ∞. El valor por defecto es p = 2, que corresponde a la longitud euclidiana o la magnitud del vector:

v = [2 0 -1];
[norm(v,1) norm(v) norm(v,inf)]
ans =

    3.0000    2.2361    2.0000

La norma-p de una matriz A,

Ap=maxxAxpxp,

se puede calcular para p = 1, 2, e ∞ usando norm(A,p). Nuevamente, el valor predeterminado es p = 2:

A = pascal(3);
[norm(A,1) norm(A) norm(A,inf)]
ans =

   10.0000    7.8730   10.0000

Cuando desee calcular la norma de cada fila o columna de una matriz, puede usar vecnorm:

vecnorm(A)
ans =

    1.7321    3.7417    6.7823

Uso de multiprocesos con funciones de álgebra lineal

MATLAB admite multiprocesos computacionales para varias funciones de álgebra lineal y funciones numéricas que actúan elemento por elemento. Estas funciones se ejecutan automáticamente en varios procesos. Para que una función o expresión se ejecute más rápido en varias CPU, se deben cumplir una serie de condiciones:

  1. La función debe realizar operaciones que se dividan fácilmente en secciones ejecutadas de manera simultánea. Estas secciones se deben poder ejecutar con poca comunicación entre los procesos. Además, estas deben requerir pocas operaciones secuenciales.

  2. El tamaño de los datos debe ser lo suficientemente grande para que cualquier ventaja propia de la ejecución simultánea justifique el tiempo que se requiere para dividir los datos y administrar procesos de ejecución separados. Por ejemplo, la mayoría de las funciones se acelera solo cuando el arreglo contiene varios miles de elementos o más.

  3. La operación no debe estar limitada por la memoria; el tiempo de procesamiento no debe estar dominado por el tiempo de acceso de la memoria. Como regla general, las funciones complicadas se aceleran más que las funciones simples.

Los operadores de multiplicación de matrices (X*Y) y potencia de matrices (X^p) muestran un importante aumento de la velocidad en arreglos de doble precisión de gran tamaño (del orden de diez mil elementos). Las funciones de análisis de matrices det, rcond, hess y expm también muestran un incremento significativo de la velocidad en arreglos de doble precisión de gran tamaño.

Temas relacionados

Sistemas de ecuaciones lineales

Consideraciones computacionales

Uno de los problemas más importantes del cálculo técnico es la solución de sistemas de ecuaciones lineales simultáneas.

En notación matricial, el problema general se plasma de la siguiente manera: Dadas dos matrices A y b, ¿existe una matriz única x de manera tal que Axb o xAb?

Resulta didáctico considerar un ejemplo de 1 por 1. Por ejemplo, la ecuación

7x = 21

¿tiene una solución única?

Naturalmente, la respuesta es sí. La solución única de la ecuación es x = 3. La solución se obtiene fácilmente a través de la división:

x = 21/7 = 3.

La solución no se obtiene comúnmente mediante el cálculo del inverso de 7, es decir 7–1= 0.142857..., y luego multiplicando 7–1 por 21. Este proceso sería más arduo, y si 7–1 se representa como un número finito de dígitos, también sería menos preciso. Pueden aplicarse consideraciones muy semejantes a los conjuntos de ecuaciones lineales con más de una incógnita. MATLAB resuelve dichas ecuaciones sin calcular inversas de matrices.

Aunque no es una notación matemática estándar, MATLAB utiliza la terminología de la división, común en el caso de escalares, para describir la solución de un sistema general de ecuaciones simultáneas. Los dos símbolos de división, barra diagonal, /, y barra invertida, \, corresponden a las dos funciones mrdivide y mldivide de MATLAB. Estos operadores se utilizan para las dos situaciones en las que la matriz desconocida aparece a la izquierda o a la derecha de la matriz de coeficientes:

x = b/A

Indica la solución a la ecuación de matrices xA = b, que se obtiene usando mrdivide.

x = A\b

Indica la solución a la ecuación de matrices Ax = b, que se obtiene usando mldivide.

Considere “dividir” ambos lados de la ecuación, Ax = b o xA = b por A. La matriz de coeficientes A siempre está en el “denominador”.

Las condiciones de compatibilidad de dimensiones para x = A\b requieren que las matrices A y b tengan la misma cantidad de filas. Entonces, la solución x tiene la misma cantidad de columnas que b y su dimensión de filas es igual a la dimensión de columnas de A. Para x = b/A, los roles de las filas y las columnas se intercambian.

En la práctica, las ecuaciones lineales del tipo Ax = b se dan con mayor frecuencia que aquellas del tipo xA = b. En consecuencia, la barra invertida se usa con mucha mayor frecuencia que la barra diagonal. El resto de esta sección se concentra en el operador barra invertida; las propiedades correspondientes del operador barra diagonal se pueden inferir de la identidad:

(b/A)' = (A'\b').

No es necesario que la matriz de coeficientes A sea cuadrada. Si A tiene un tamaño de m por n, entonces, existen tres casos:

m = n

Sistema cuadrado. Busque una solución exacta.

m > n

Sistema sobredeterminado, con más ecuaciones que incógnitas. Encuentre una solución de mínimos cuadrados.

m < n

Sistema subdeterminado, con menos ecuaciones que incógnitas. Encuentre una solución básica con un máximo de m componentes distintos de cero.

El algoritmo de mldivide.  El operador mldivide emplea diferentes mecanismos de solución para tratar distintos tipos de matrices de coeficientes. Los diversos casos se diagnostican de forma automática mediante el análisis de la matriz de coeficientes. Para obtener más información, consulte la sección “Algoritmos” de la página de referencia de mldivide.

Solución general

La solución general para un sistema de ecuaciones lineales Axb describe todas las soluciones posibles. Se puede encontrar la solución general mediante:

  1. La solución del sistema homogéneo correspondiente Ax = 0. Para hacerlo, use el comando null, escribiendo null(A). Esto devuelve una base para el espacio de solución de Ax = 0. Cualquier solución es una combinación lineal de vectores de base.

  2. La búsqueda de una solución particular para el sistema no homogéneo Ax =b.

Luego, se puede escribir cualquier solución de Axb como la suma de la solución particular para Ax =b, del paso 2, más una combinación lineal de los vectores de base del paso 1.

En el resto de esta sección, se describe cómo usar MATLAB para buscar una solución particular para Ax =b, como en el paso 2.

Sistemas cuadrados

La situación más común incluye una matriz de coeficientes cuadrada A y un único vector columna b a la derecha.

Matriz de coeficientes no singular.  Si la matriz A no es singular, entonces la solución, x = A\b, tiene el mismo tamaño que b. Por ejemplo:

A = pascal(3);
u = [3; 1; 4];
x = A\u

x =
      10
     -12
       5

Puede confirmarse que A*x es exactamente igual a u.

Si A y b son cuadradas y tienen el mismo tamaño, x= A\b también tiene ese tamaño:

b = magic(3);
X = A\b

X =
      19    -3    -1
     -17     4    13
       6     0    -6

Puede confirmarse que A*x es exactamente igual a b.

Estos dos ejemplos tienen soluciones exactas y de números enteros. Esto se debe a que se escogió como matriz de coeficientes pascal(3), que es una matriz de rango completo (no singular).

Matriz de coeficientes singular.  Una matriz cuadrada A es singular si no tiene columnas linealmente independientes. Si A es singular, la solución para Ax = b o bien no existe o no es única. El operador barra invertida, A\b, emite una advertencia si A es casi singular o si detecta una singularidad exacta.

Si A es singular y Ax = b tiene solución, se puede encontrar una solución particular que no es única. Para ello, escriba

P = pinv(A)*b

pinv(A) es una pseudoinversa de A. Si Ax = b no tiene solución exacta, entonces pinv(A) devuelve una solución de mínimos cuadrados.

Por ejemplo:

A = [ 1     3     7
     -1     4     4
      1    10    18 ]

es singular, como puede verificar si escribe

rank(A)

ans =

     2

Debido a que A no es de rango completo, tiene algunos valores singulares iguales a cero.

Soluciones exactas Para b =[5;2;12], la ecuación Ax = b tiene una solución exacta, dada por:

pinv(A)*b

ans =
    0.3850
   -0.1103
    0.7066

Para comprobar que pinv(A)*b es una solución exacta, escriba

A*pinv(A)*b

ans =
    5.0000
    2.0000
   12.0000

Soluciones de mínimos cuadrados. Sin embargo, si b = [3;6;0], Ax = b no tiene solución exacta. En este caso, pinv(A)*b devuelve una solución de mínimos cuadrados. Si escribe

A*pinv(A)*b

ans =
   -1.0000
    4.0000
    2.0000

no recupera el vector original b.

Puede determinar si Ax =b tiene una solución exacta buscando la forma escalonada reducida por filas de la matriz aumentada [A b]. Para hacerlo en este ejemplo, introduzca

rref([A b])
ans =
    1.0000         0    2.2857         0
         0    1.0000    1.5714         0
         0         0         0    1.0000

Dado que la fila inferior contiene solo ceros, con excepción de la última entrada, la ecuación no tiene solución. En este caso, pinv(A) devuelve una solución de mínimos cuadrados.

Sistemas sobredeterminados

Este ejemplo muestra cómo con frecuencia se encuentran sistemas sobredeterminados en diversos tipos de ajustes de curvas para datos experimentales.

Una cantidad, y, se mide en varios valores de tiempo diferentes, t, para producir las siguientes observaciones. Para introducir los datos y verlos en una tabla, utilice las siguientes instrucciones.

t = [0 .3 .8 1.1 1.6 2.3]';
y = [.82 .72 .63 .60 .55 .50]';
B = table(t,y)
B=6×2 table
     t      y  
    ___    ____

      0    0.82
    0.3    0.72
    0.8    0.63
    1.1     0.6
    1.6    0.55
    2.3     0.5

Intente modelar los datos con una función exponencial decreciente

y(t)=c1+c2e-t.

La ecuación anterior indica que el vector y se debe aproximar mediante una combinación lineal de otros dos vectores. Uno es un vector constante que contiene solo unos, y el otro es el vector con componentes exp(-t). Los coeficientes desconocidos, c1 y c2, se pueden calcular haciendo un ajuste de mínimos cuadrados que minimice la suma de los cuadrados de las desviaciones de los datos respecto del modelo. Existen seis ecuaciones en dos incógnitas, representadas mediante una matriz de 6 por 2.

E = [ones(size(t)) exp(-t)]
E = 6×2

    1.0000    1.0000
    1.0000    0.7408
    1.0000    0.4493
    1.0000    0.3329
    1.0000    0.2019
    1.0000    0.1003

Use el operador barra invertida para obtener la solución de mínimos cuadrados.

c = E\y
c = 2×1

    0.4760
    0.3413

En otras palabras, el ajuste de mínimos cuadrados para los datos es el siguiente:

y(t)=0.4760+0.3413e-t.

Las siguientes instrucciones evalúan el modelo en incrementos uniformemente espaciados en t, y luego representan gráficamente el resultado junto con los datos originales:

T = (0:0.1:2.5)';
Y = [ones(size(T)) exp(-T)]*c;
plot(T,Y,'-',t,y,'o')

E*c no es exactamente igual a y, pero la diferencia bien podría ser menor que los errores de medición en los datos originales.

Una matriz rectangular A es de rango deficiente si no tiene columnas linealmente independientes. Si A es de rango deficiente, entonces, no hay una solución única de mínimos cuadrados para AX = B. A\B emite una advertencia si A es de rango deficiente y produce una solución de mínimos cuadrados. Puede utilizar lsqminnorm para encontrar la solución X que tiene la norma mínima entre todas las soluciones.

Sistemas subdeterminados

En este ejemplo se demuestra cómo la solución para los sistemas subdeterminados no es única. Los sistemas lineales subdeterminados incluyen más incógnitas que ecuaciones. En MATLAB, la operación de división matricial izquierda encuentra una solución de mínimos cuadrados básica que, como máximo, tiene m componentes distintos de cero para una matriz de coeficientes de m por n.

El siguiente es un breve ejemplo al azar:

R = [6 8 7 3; 3 5 4 1]
rng(0);
b = randi(8,2,1)
R =

       6              8              7              3       
       3              5              4              1       


b =

       7       
       8      

El sistema lineal Rp = b incluye dos ecuaciones en cuatro incógnitas. Como la matriz de coeficientes contiene enteros pequeños, es apropiado usar el comando format para visualizar la solución en formato racional. La solución particular se obtiene de la siguiente manera:

format rat
p = R\b
p =

       0       
      17/7     
       0       
     -29/7    

Uno de los componentes distintos de cero es p(2), debido a que R(:,2) es la columna de R con la norma de mayor tamaño. El otro componente que no es cero es p(4), ya que R(:,4) domina tras la eliminación de R(:,2).

La solución general completa al sistema subdeterminado se puede caracterizar mediante la adición de p a una combinación lineal arbitraria de los vectores del espacio nulo, que se pueden encontrar usando la función null con una opción que solicita una base racional.

Z = null(R,'r')
Z =

      -1/2           -7/6     
      -1/2            1/2     
       1              0       
       0              1       

Se puede confirmar que R*Z es igual a cero y que el residuo R*x - b resulta pequeño para cualquier vector x, donde

x = p + Z*q

Dado que las columnas de Z son los vectores del espacio nulo, el producto Z*q es una combinación lineal de dichos vectores:

Zq=(x1x2)(uw)=ux1+wx2.

Para ilustrarlo, escoja un q arbitrario y construya x.

q = [-2; 1];
x = p + Z*q;

Calcule la norma del residuo.

format short
norm(R*x - b)
ans =

   2.6645e-15

Cuando hay infinitas soluciones disponibles, la solución con la norma mínima resulta de especial interés. Puede utilizar lsqminnorm para calcular la solución de mínimos cuadrados de norma mínima. Esta solución tiene el menor valor posible para norm(p).

p = lsqminnorm(R,b)
p =

    -207/137   
     365/137   
      79/137   
    -424/137  

Resolución de varios lados derechos

Algunos problemas se relacionan con la resolución de sistemas lineales que tienen la misma matriz de coeficientes A, pero diferentes lados derechos b. Cuando los diferentes valores de b están disponibles al mismo tiempo, puede construir b como matriz con varias columnas y resolver a la vez todos los sistemas de ecuaciones mediante un único comando de barra invertida: X = A\[b1 b2 b3 …].

No obstante, a veces, todos los valores de b no están disponibles al mismo tiempo, lo que significa que tendrá que resolver varios sistemas de ecuaciones de forma consecutiva. Cuando resuelve uno de estos sistemas de ecuaciones mediante la barra diagonal (/) o la barra invertida (\), el operador factoriza la matriz de coeficientes A y utiliza esta descomposición de matrices para calcular la solución. Sin embargo, a partir de entonces, cada vez que resuelva un sistema de ecuaciones similar con una b diferente, el operador calculará la misma descomposición de A, lo cual es un cálculo redundante.

La solución para este problema es realizar un cálculo previo con la descomposición de A y, a continuación, reutilizar los factores para resolver los diferentes valores de b. Sin embargo, en la práctica, puede que sea difícil realizar el cálculo previo de la descomposición de esta manera, ya que se debe saber qué descomposición calcular (LU, LDL, Cholesky, etc.) y cómo multiplicar los factores para resolver el problema. Por ejemplo, con la descomposición LU, se deben resolver dos sistemas lineales para resolver el sistema original Ax = b:

[L,U] = lu(A);
x = U \ (L \ b);

En cambio, para resolver sistemas lineales con varios lados derechos consecutivos, lo más recomendado es usar objetos decomposition. Estos objetos le permiten aprovechar las ventajas de rendimiento que ofrece el cálculo previo de la descomposición de matrices, pero no es necesario saber cómo usar los factores de matrices. Puede reemplazar la descomposición de LU anterior de la siguiente manera:

dA = decomposition(A,'lu');
x = dA\b;

Si no sabe con certeza qué descomposición utilizar, decomposition(A) escoge el tipo correcto en función de las propiedades de A, algo similar a lo que hace la barra invertida.

Aquí le mostramos una prueba simple sobre las posibles ventajas de rendimiento de este enfoque. La prueba resuelve el mismo sistema lineal disperso 100 veces usando la barra invertida (\) y la función decomposition.

n = 1e3;
A = sprand(n,n,0.2) + speye(n);
b = ones(n,1);

% Backslash solution
tic
for k = 1:100
    x = A\b;
end
toc
Elapsed time is 9.006156 seconds.
% decomposition solution
tic
dA = decomposition(A);
for k = 1:100
    x = dA\b;
end
toc
Elapsed time is 0.374347 seconds.

Para este problema, la solución decomposition es mucho más rápida que la ofrecida por el uso exclusivo de la barra invertida. Sin embargo, la sintaxis sigue siendo simple.

Métodos iterativos

Si la matriz de coeficientes A es grande y dispersa, generalmente, los métodos de factorización no son eficientes. Los métodos iterativos generan una serie de soluciones aproximadas. MATLAB ofrece varios métodos iterativos para manipular matrices de entrada grandes y dispersas.

FunciónDescripción
pcg

Método del gradiente conjugado precondicionado Este método es adecuado para la matriz de coeficientes definida positiva hermítica A.

bicg

Método del gradiente biconjugado

bicgstab

Método del gradiente biconjugado estabilizado

bicgstabl

Método BiCGStab(l), variante del método del gradiente biconjugado estabilizado

cgs

Método del gradiente conjugado cuadrado

gmres

Método del residuo mínimo generalizado

lsqr

Método LSQR

minres

Método del residuo mínimo Este método es adecuado para la matriz de coeficientes hermítica A.

qmr

Método del residuo cuasi-mínimo

symmlq

Método LQ simétrico

tfqmr

Método QMR sin trasposición

Multiproceso computacional

MATLAB admite multiprocesos computacionales para varias funciones de álgebra lineal y funciones numéricas que actúan elemento por elemento. Estas funciones se ejecutan automáticamente en varios procesos. Para que una función o expresión se ejecute más rápido en varias CPU, se deben cumplir una serie de condiciones:

  1. La función debe realizar operaciones que se dividan fácilmente en secciones ejecutadas de manera simultánea. Estas secciones se deben poder ejecutar con poca comunicación entre los procesos. Además, estas deben requerir pocas operaciones secuenciales.

  2. El tamaño de los datos debe ser lo suficientemente grande para que cualquier ventaja propia de la ejecución simultánea justifique el tiempo que se requiere para dividir los datos y administrar procesos de ejecución separados. Por ejemplo, la mayoría de las funciones se acelera solo cuando el arreglo contiene varios miles de elementos o más.

  3. La operación no debe estar limitada por la memoria; el tiempo de procesamiento no debe estar dominado por el tiempo de acceso de la memoria. Como regla general, las funciones complicadas se aceleran más que las funciones simples.

Las funciones inv, lscov, linsolve y mldivide muestran un significativo aumento de la velocidad en arreglos de doble precisión de gran tamaño (del orden de diez mil elementos o más) cuando está activado el multiproceso.

Consulte también

| | | |

Temas relacionados

Factorizaciones

Introducción

Las tres factorizaciones de matrices que se analizan en esta sección utilizan matrices triangulares, donde todos los elementos ubicados por encima o por debajo de la diagonal son cero. Los sistemas de ecuaciones lineales que incluyen matrices triangulares se resuelven fácil y rápidamente a través de la sustitución hacia adelante o hacia atrás.

Factorización de Cholesky

La factorización de Cholesky expresa una matriz simétrica como el producto de una matriz triangular y su matriz traspuesta.

A = RR,

donde R es una matriz triangular superior.

No todas las matrices simétricas se pueden factorizar de esta manera. Las matrices que tienen dicha factorización se califican como definidas positivas. Esto implica que todos los elementos diagonales de A son positivos y que los elementos fuera de la diagonal “no son demasiado grandes”. Las matrices de Pascal proporcionan un ejemplo interesante. A lo largo de este capítulo, la matriz de ejemplo A ha sido la matriz de Pascal de 3 por 3. Por un momento, considere la matriz de 6 por 6:

A = pascal(6)

A =
       1     1     1     1     1     1
       1     2     3     4     5     6
       1     3     6    10    15    21
       1     4    10    20    35    56
       1     5    15    35    70   126
       1     6    21    56   126   252

Los elementos de A son coeficientes binomiales. Cada elemento es la suma de sus vecinos del norte y el oeste. La factorización de Cholesky es la siguiente:

R = chol(A)

R =
     1     1     1     1     1     1
     0     1     2     3     4     5
     0     0     1     3     6    10
     0     0     0     1     4    10
     0     0     0     0     1     5
     0     0     0     0     0     1

Los elementos son, nuevamente, coeficientes binomiales. El hecho de que R'*R es igual a A demuestra una identidad que incluye las sumas de productos de coeficientes binomiales.

Nota

La factorización de Cholesky también se aplica a matrices complejas. Cualquier matriz compleja que tiene una factorización de Cholesky satisface la condición

A′ = A

y se califica como definida positiva hermítica.

La factorización de Cholesky permite que el sistema lineal

Ax = b

sea reemplazado por

RRx = b.

Debido a que el operador barra invertida reconoce sistemas triangulares, este sistema puede resolverse fácilmente en el entorno de MATLAB con la siguiente operación:

x = R\(R'\b)

Si A es n por n, la complejidad computacional de chol(A) es O(n3), pero la complejidad de las soluciones de barra invertida subsiguientes es solo O(n2).

Factorización LU

La factorización LU, o eliminación gaussiana, expresa cualquier matriz cuadrada A como el producto de una permutación de una matriz triangular inferior y una matriz triangular superior:

A = LU,

donde L es una permutación de una matriz triangular inferior con unos en su diagonal y U es una matriz triangular superior.

Las permutaciones son necesarias por motivos tanto teóricos como computacionales. La matriz

[0110]

no se puede expresar como el producto de matrices triangulares sin intercambiar sus dos filas. Aunque la matriz

[ε110]

se puede expresar como el producto de matrices triangulares, cuando ε es pequeño, los elementos de los factores son grandes y amplían los errores. Por lo tanto, aunque las permutaciones no son estrictamente necesarias, resultan convenientes. El pivoteo parcial asegura que los elementos de L están acotados por uno en magnitud y que los elementos de U no son mucho más grandes que los de A.

Por ejemplo:

[L,U] = lu(B)

L =
    1.0000         0         0
    0.3750    0.5441    1.0000
    0.5000    1.0000         0

U =
    8.0000    1.0000    6.0000
         0    8.5000   -1.0000
         0         0    5.2941

La factorización LU de A permite que el sistema lineal

A*x = b

se resuelva rápidamente de la siguiente manera:

x = U\(L\b)

Los determinantes y las inversas se calculan a partir de la factorización LU con la fórmula

det(A) = det(L)*det(U)

y la fórmula

inv(A) = inv(U)*inv(L)

También es posible calcular los determinantes con det(A) = prod(diag(U)), aunque los signos de los determinantes podrían invertirse.

Factorización QR

Una matriz ortogonal, o matriz con columnas ortogonales, es una matriz real cuyas columnas tienen longitud uno y son perpendiculares entre sí. Si Q es ortogonal, entonces

QQ = 1.

Las matrices ortogonales más simples son las rotaciones de coordenadas bidimensionales:

[cos(θ)sin(θ)sin(θ)cos(θ)].

En el caso de las matrices complejas, el término correspondiente es unitaria. Las matrices ortogonales y unitarias son convenientes para el cálculo numérico porque mantienen la longitud, conservan los ángulos y no amplían los errores.

La factorización ortogonal, o QR, expresa cualquier matriz rectangular como el producto de una matriz ortogonal o unitaria y una matriz triangular superior.También podría incluir una permutación de columnas:

A = QR

o

AP = QR,

donde Q es ortogonal o unitaria, R es triangular superior y P es una permutación.

Existen cuatro variantes de la factorización QR: de tamaño total o parcial, y con permutación de columnas o sin ella.

Los sistemas lineales sobredeterminados incluyen una matriz rectangular con más filas que columnas, es decir, m por n con m > n. La factorización QR de tamaño total produce una matriz cuadrada ortogonal Q de m por m y una matriz rectangular triangular superior R de m por n:

C=gallery('uniformdata',[5 4], 0);
[Q,R] = qr(C)

Q =

    0.6191    0.1406   -0.1899   -0.5058    0.5522
    0.1506    0.4084    0.5034    0.5974    0.4475
    0.3954   -0.5564    0.6869   -0.1478   -0.2008
    0.3167    0.6676    0.1351   -0.1729   -0.6370
    0.5808   -0.2410   -0.4695    0.5792   -0.2207


R =

    1.5346    1.0663    1.2010    1.4036
         0    0.7245    0.3474   -0.0126
         0         0    0.9320    0.6596
         0         0         0    0.6648
         0         0         0         0

En muchos casos, las últimas m – n columnas de Q no son necesarias, ya que se multiplican por los ceros de la parte inferior de R. Por lo tanto, la factorización QR de tamaño parcial produce una matriz rectangular Q de m por n con columnas ortonormales y una matriz cuadrada triangular superior R de n por n. Para el ejemplo de 5 por 4, el ahorro no es tan significativo, pero para matrices más grandes y altamente rectangulares, el ahorro tanto de tiempo como de memoria puede ser bastante importante:

[Q,R] = qr(C,0)	
Q =

    0.6191    0.1406   -0.1899   -0.5058
    0.1506    0.4084    0.5034    0.5974
    0.3954   -0.5564    0.6869   -0.1478
    0.3167    0.6676    0.1351   -0.1729
    0.5808   -0.2410   -0.4695    0.5792


R =

    1.5346    1.0663    1.2010    1.4036
         0    0.7245    0.3474   -0.0126
         0         0    0.9320    0.6596
         0         0         0    0.6648

En contraste con la factorización LU, la factorización QR no requiere pivoteo ni permutaciones. Sin embargo, es de utilidad una permutación de columnas opcional, activada por la presencia de un tercer argumento de salida, para poder detectar singularidad o deficiencia de rango. En cada paso de la factorización, se usa como base la columna con la norma más grande de la matriz restante sin factorizar. Esto asegura que los elementos diagonales de R aparezcan en orden descendente y que cualquier dependencia lineal entre las columnas se pueda revelar, casi con toda certeza, al observar estos elementos. Para el pequeño ejemplo que se muestra aquí, la segunda columna de C tiene una norma más grande que la primera, por lo que las dos columnas se intercambian:

[Q,R,P] = qr(C)

Q =
   -0.3522    0.8398   -0.4131
   -0.7044   -0.5285   -0.4739
   -0.6163    0.1241    0.7777

R =
  -11.3578   -8.2762
         0    7.2460
         0         0

P =
     0     1
     1     0

Cuando se combinan el tamaño parcial y las permutaciones de columna, el tercer argumento de salida es un vector de permutación, y no una matriz de permutación:

[Q,R,p] = qr(C,0)

Q =
   -0.3522    0.8398
   -0.7044   -0.5285
   -0.6163    0.1241

R =
  -11.3578   -8.2762
         0    7.2460


p =
     2     1

La factorización QR transforma un sistema lineal sobredeterminado en un sistema triangular equivalente. La expresión

norm(A*x - b)

es igual a

norm(Q*R*x - b)

La multiplicación por matrices ortogonales conserva la norma euclidiana, de manera que esta expresión también es igual a

norm(R*x - y)

donde y = Q'*b. Dado que las últimas filas m-n de R son cero, esta expresión se divide en dos partes:

norm(R(1:n,1:n)*x - y(1:n))

y

norm(y(n+1:m))

Cuando A tiene rango completo, es posible despejar la x de manera que la primera de estas expresiones sea igual a cero. Entonces, la segunda expresión devuelve la norma del residuo. Cuando A no tiene rango completo, la estructura triangular de R permite hallar una solución básica para el problema de mínimos cuadrados.

Uso de multiprocesos para factorización

El software de MATLAB admite multiprocesos computacionales para varias funciones de álgebra lineal y funciones numéricas que actúan elemento por elemento. Estas funciones se ejecutan automáticamente en varios procesos. Para que una función o expresión se ejecute más rápido en varias CPU, se deben cumplir una serie de condiciones:

  1. La función debe realizar operaciones que se dividan fácilmente en secciones ejecutadas de manera simultánea. Estas secciones se deben poder ejecutar con poca comunicación entre los procesos. Además, estas deben requerir pocas operaciones secuenciales.

  2. El tamaño de los datos debe ser lo suficientemente grande para que cualquier ventaja propia de la ejecución simultánea justifique el tiempo que se requiere para dividir los datos y administrar procesos de ejecución separados. Por ejemplo, la mayoría de las funciones se acelera solo cuando el arreglo contiene varios miles de elementos o más.

  3. La operación no debe estar limitada por la memoria; el tiempo de procesamiento no debe estar dominado por el tiempo de acceso de la memoria. Como regla general, las funciones complicadas se aceleran más que las funciones simples.

Las funciones lu y qr muestran un importante incremento de la velocidad en arreglos de precisión doble de gran tamaño (del orden de diez mil elementos).

Consulte también

| |

Temas relacionados

Potencias y exponenciales

Potencias enteras positivas

Si A es una matriz cuadrada y p es un entero positivo, A^p efectivamente multiplica a A por sí misma p-1 veces. Por ejemplo:

A = [1 1 1;1 2 3;1 3 6]

A =

     1     1     1
     1     2     3
     1     3     6

X = A^2

X =
     3     6    10
     6    14    25
    10    25    46

Potencias inversas y racionales

Si A es cuadrada y no singular, A^(-p) efectivamente multiplica a inv(A) por sí misma p-1 veces:

Y = A^(-3)

Y =

  145.0000 -207.0000   81.0000
 -207.0000  298.0000 -117.0000
   81.0000 -117.0000   46.0000

También se permiten las potencias racionales, como A^(2/3), y los resultados dependen de la distribución de los valores propios de la matriz.

Potencias elemento por elemento

El operador .^ produce potencias elemento por elemento. Por ejemplo:

X = A.^2

A =
     1     1     1
     1     4     9
     1     9    36

Exponenciales

La función sqrtm(A) calcula A^(1/2) mediante un algoritmo más exacto. La m en sqrtm distingue esta función de sqrt(A), la cual, al igual que A.^(1/2), opera elemento por elemento.

Un sistema de ecuaciones diferenciales lineales de coeficientes constantes se puede escribir como

dx/dt=Ax,

donde x=x(t) es un vector de funciones de t y A es una matriz independiente de t. La solución se puede expresar en términos de la exponencial de matriz

x(t)=etAx(0).

La función expm(A) calcula la exponencial de una matriz. Un ejemplo lo proporcionan la matriz de coeficientes de 3 por 3,

A = [0 -6 -1; 6 2 -16; -5 20 -10]
A = 3×3

     0    -6    -1
     6     2   -16
    -5    20   -10

y la condición inicial, x(0).

x0 = [1 1 1]'
x0 = 3×1

     1
     1
     1

La exponencial de una matriz se utiliza para calcular la solución, x(t), a la ecuación diferencial en 101 puntos en el intervalo 0t1.

X = [];
for t = 0:.01:1
   X = [X expm(t*A)*x0]; 
end

Un diagrama de plano de fase tridimensional muestra la solución en trayectoria espiral hacia el origen. Este comportamiento está relacionado con los valores propios de la matriz de coeficientes.

plot3(X(1,:),X(2,:),X(3,:),'-o')

Consulte también

| |

Temas relacionados

Valores propios

Descomposición en valores propios

Un valor propio y un vector propio de una matriz cuadrada A son, respectivamente, un escalar λ y un vector distinto de cero υ que satisfacen

= λυ.

Con los valores propios en la diagonal de una matriz diagonal Λ y los correspondientes vectores propios formando las columnas de una matriz V, se tiene que

AV = .

Si V no es singular, esta expresión se convierte en la descomposición en valores propios

A = VΛV–1.

Un buen ejemplo es la matriz de coeficientes de esta ecuación diferencial ordinaria:

A =
     0    -6    -1
     6     2   -16
    -5    20   -10

La instrucción

lambda = eig(A)

produce un vector columna que contiene los valores propios. Para esta matriz, los valores propios son complejos:

lambda =
     -3.0710         
     -2.4645+17.6008i
     -2.4645-17.6008i

La parte real de cada uno de los valores propios es negativa, y por eso eλt tiende a cero a medida que aumenta t. La parte imaginaria no nula de dos de los valores propios, ±ω, aporta el componente oscilatorio, sin(ωt), a la solución de la ecuación diferencial.

Con dos argumentos de salida, eig calcula los vectores propios y almacena los valores propios en una matriz diagonal:

[V,D] = eig(A)

V =
  -0.8326         0.2003 - 0.1394i   0.2003 + 0.1394i
  -0.3553        -0.2110 - 0.6447i  -0.2110 + 0.6447i
  -0.4248        -0.6930            -0.6930          

D =
  -3.0710                 0                 0         
        0           -2.4645+17.6008i        0         
        0                 0           -2.4645-17.6008i

El primer vector propio es real y los otros dos vectores son conjugados complejos entre sí. Los tres vectores están normalizados de manera que su longitud euclidiana, norm(v,2), es igual a uno.

La matriz V*D*inv(V), que se puede escribir de forma más concisa como V*D/V, coincide con A dentro del error de redondeo. Por su parte, inv(V)*A*V, o V\A*V, coincide con D dentro del error de redondeo.

Valores propios múltiples

Algunas matrices no tienen descomposición en valores propios. Estas matrices no se pueden diagonalizar. Por ejemplo:

A = [ 1    -2    1 
      0     1    4 
      0     0    3 ]

Para esta matriz

[V,D] = eig(A)

produce

V =

    1.0000    1.0000   -0.5571
         0    0.0000    0.7428
         0         0    0.3714


D =

     1     0     0
     0     1     0
     0     0     3

Hay un valor propio doble en λ = 1. La primera y la segunda columna de V son iguales. Para esta matriz, no existe un conjunto completo de vectores propios linealmente independientes.

Descomposición de Schur

Muchos cálculos matriciales avanzados no requieren descomposiciones en valores propios. En su lugar, están basados en la descomposición de Schur

A = USU ′ ,

donde U es una matriz ortogonal y S es una matriz triangular superior por bloques, con bloques de 1 por 1 y de 2 por 2 en la diagonal. Los valores propios se muestran mediante los elementos y bloques diagonales de S, mientras que las columnas de U proporcionan una base ortogonal con propiedades numéricas mucho mejores que las de un conjunto de vectores propios.

Por ejemplo, compare el valor propio y las descomposiciones Schur de esta matriz defectuosa:

A = [ 6    12    19 
     -9   -20   -33 
      4     9    15 ];

[V,D] = eig(A)

V =

  -0.4741 + 0.0000i  -0.4082 - 0.0000i  -0.4082 + 0.0000i
   0.8127 + 0.0000i   0.8165 + 0.0000i   0.8165 + 0.0000i
  -0.3386 + 0.0000i  -0.4082 + 0.0000i  -0.4082 - 0.0000i


D =

  -1.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   1.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   1.0000 - 0.0000i

[U,S] = schur(A)

U =

   -0.4741    0.6648    0.5774
    0.8127    0.0782    0.5774
   -0.3386   -0.7430    0.5774


S =

   -1.0000   20.7846  -44.6948
         0    1.0000   -0.6096
         0    0.0000    1.0000

La matriz A es defectuosa, ya que no tiene un conjunto completo de vectores propios linealmente independientes (la segunda y tercera columna de V son iguales). Debido a que no todas las columnas de V son linealmente independientes, tiene un número de condición grande de aproximadamente ~1e8. Sin embargo, schur puede calcular tres vectores de base diferentes en U. En vista de que U es ortogonal, cond(U) = 1.

La matriz S tiene el valor propio real como la primera entrada en la diagonal y el valor propio repetido representado por el bloque de 2 por 2 de la parte inferior derecha. Los valores propios del bloque de 2 por 2 también son valores propios de A:

eig(S(2:3,2:3))

ans =

   1.0000 + 0.0000i
   1.0000 - 0.0000i

Consulte también

|

Temas relacionados

Valores singulares

Un valor singular y los correspondientes vectores singulares de una matriz rectangular A son, respectivamente, un escalar σ y un par de vectores u y v que satisfacen

Av=σuAHu=σv,

donde AH es la traspuesta hermítica de A. Los vectores singulares u y v comúnmente se normalizan a 1. Igualmente, si u y v son vectores singulares de A, entonces -u y -v también son vectores singulares de A.

Los valores singulares σ siempre son reales y no negativos, aunque A sea compleja. Con los valores singulares en la diagonal de una matriz diagonal Σ y los correspondientes vectores singulares formando las columnas de dos matrices ortogonales U y V, se obtienen las ecuaciones

AV=UΣAHU=VΣ.

Como U y V son matrices unitarias, si se multiplica la primera ecuación por VH a la derecha, se produce la ecuación de descomposición en valores singulares

A=UΣVH.

La descomposición en valores singulares completa de una matriz de m por n incluye una U de m por m, una Σ de m por n y una V de n por n. En otras palabras, U y V son cuadradas y Σ tiene el mismo tamaño que A. Si A tiene muchas más filas que columnas (m > n), entonces la matriz U de m por m resultante es grande. Sin embargo, la mayoría de las columnas de U se multiplica por ceros en Σ. En esta situación, la descomposición de tamaño parcial ahorra tanto tiempo como almacenamiento, ya que produce una U de m por n, una Σ de n por n y la misma V:

La descomposición en valores propios es la herramienta apropiada para analizar una matriz cuando esta representa una aplicación de un espacio vectorial en sí mismo, como en el caso de una ecuación diferencial ordinaria. Sin embargo, la descomposición en valores singulares es la herramienta adecuada para analizar una aplicación de un espacio vectorial en otro espacio vectorial, posiblemente con una dimensión diferente. La mayoría de los sistemas de ecuaciones lineales simultáneas encaja en esta segunda categoría.

Si A es cuadrada, simétrica y definida positiva, entonces sus descomposiciones en valores propios y valores singulares son iguales. Sin embargo, a medida que A se desvía de la simetría y definición positiva, la diferencia entre las dos descomposiciones aumenta. En particular, la descomposición en valores singulares de una matriz real siempre es real, pero la descomposición de valores propios de una matriz real no simétrica puede ser compleja.

Para la matriz ejemplo

A =
     9     4
     6     8
     2     7

la descomposición en valores singulares completa es

[U,S,V] = svd(A)

U =

    0.6105   -0.7174    0.3355
    0.6646    0.2336   -0.7098
    0.4308    0.6563    0.6194


S =

   14.9359         0
         0    5.1883
         0         0


V =

    0.6925   -0.7214
    0.7214    0.6925

Se puede verificar que U*S*V' es igual a A dentro del error de redondeo. Para este pequeño problema, la descomposición de tamaño parcial es solo ligeramente menor.

[U,S,V] = svd(A,0)

U =

    0.6105   -0.7174
    0.6646    0.2336
    0.4308    0.6563


S =

   14.9359         0
         0    5.1883


V =

    0.6925   -0.7214
    0.7214    0.6925

Nuevamente, U*S*V' es igual a A dentro del error de redondeo.

Cuando la matriz A es grande y dispersa, no siempre es práctico usar svd para calcular todos los vectores y valores singulares. Por ejemplo, si necesita conocer solo algunos de los valores singulares más grandes, implica mucho trabajo adicional calcular todos los valores singulares de una matriz dispersa de 5000 por 5000. En los casos en que se requiere solo un subconjunto de los vectores y valores singulares, se prefiere la función svds en lugar de svd.

Para una matriz dispersa aleatoria de 1000 por 1000 con una densidad de aproximadamente 30%,

n = 1000;
A = sprand(n,n,0.3);

los seis valores singulares más grandes son

S = svds(A)

S =

  130.2184
   16.4358
   16.4119
   16.3688
   16.3242
   16.2838

Además, los valores singulares más pequeños son

S = svds(A,6,'smallest')

S =

    0.0740
    0.0574
    0.0388
    0.0282
    0.0131
    0.0066

Para matrices más pequeñas que pueden caber en la memoria en su versión completa, full(A), el uso de svd(full(A)) puede ser todavía más rápido que svds. No obstante, para matrices verdaderamente grandes y dispersas, se hace necesario el uso de svds.

Consulte también

| |

Temas relacionados