2012-06-02 14 views
5

Tengo una matriz con 3 dimensiones Y (i, j, w). Quiero obtener un vector determinante d (w), en el cual cada número sería el determinante de la matriz Y (:,:, w).return determinant vector - Matlab

¿Hay una sintaxis elegante para él, o simplemente tengo que usar un bucle?

gracias

Respuesta

7

Bueno, primero de todo es que casi nunca realmente quieren para calcular un factor determinante, que acaba de pensar que haces. De hecho, casi nunca es algo bueno, porque los determinantes están escasamente dimensionados. Con demasiada frecuencia se usan para inferir el estado de singularidad de una matriz, lo cual es algo terrible de hacer en términos de análisis numérico.

haber declarado mi mini-perorata contra determinantes en general ...

OPCIÓN 1:

Convierte la matriz 3-D en una serie de células de matrices cuadradas, con cada plano de la matriz como una celda. mat2cell hará el truco de manera fácil y eficiente.

A continuación, use cellfun en la matriz de celdas. cellfun puede aplicar una función (@det) a cada celda, y luego devolverá un vector de determinantes. ¿Es esto increíblemente eficiente? Probablemente no sea una gran ventaja sobre la aplicación de det en un bucle, siempre que preasigne el vector de antemano cuando escriba un bucle.

OPCIÓN 2:

Si las matrices son pequeñas, por lo que dicen 2x2 o 3x3 matrices, y luego se expanden a cabo las multiplicaciones para el determinante se multiplica vector como explícitos. Creo que esto no es clara ya que estoy escribiendo, así que para un caso de 2x2, donde se 2x2xn Y:

d = Y(1,1,:).*Y(2,2,:) - Y(1,2,:).*Y(2,1,:); 

Si ves que esto forma un vector de determinantes 2x2 para cada plano de la matriz Y. El El caso 3x3 es lo suficientemente simple como para escribir también, como seis productos de términos tripartitos. No he revisado cuidadosamente el caso 3x3 a continuación, pero debe estar cerca.

d = Y(1,1,:).*Y(2,2,:).*Y(3,3,:) + ... 
    Y(2,1,:).*Y(3,2,:).*Y(1,3,:) + ... 
    Y(3,1,:).*Y(1,2,:).*Y(2,3,:) - ... 
    Y(3,1,:).*Y(2,2,:).*Y(1,3,:) - ... 
    Y(2,1,:).*Y(1,2,:).*Y(3,3,:) - ... 
    Y(1,1,:).*Y(3,2,:).*Y(2,3,:); 

Como puede ver, la OPCIÓN 2 será bastante rápida, y está vectorizada.

Editar: como respuesta a Chris, hay una diferencia SIGNIFICATIVA en el tiempo requerido. Considere el tiempo requerido para un conjunto de matrices 1e5.

p = 2; 
n = 1e5; 
Y = rand(p,p,n); 

tic, 
d0 = squeeze(Y(1,1,:).*Y(2,2,:) - Y(2,1,:).*Y(1,2,:)); 
toc 

Elapsed time is 0.002141 seconds. 

tic, 
X = squeeze(mat2cell(Y,p,p,ones(1,n))); 
d1= cellfun(@det,X); 
toc 

Elapsed time is 12.041883 seconds. 

Las dos llamadas devuelven los mismos valores dentro de la basura flotante.

std(d0-d1) 
ans = 
    3.8312e-17 

Un bucle no sería mejor, de hecho, seguramente peor. Entonces debía escribir un fragmento de código que se enfrentaría con la tarea de generar determinantes para MUCHAS de tales matrices en una matriz, sería un caso especial el código para las matrices 2x2 y 3x3. Incluso podría escribirlo para matrices 4x4. Sí, es un desastre escribir, pero hay una gran diferencia en el tiempo requerido.

Una razón es que el det de MATLAB utiliza una llamada a LU, factorizando la matriz. En teoría, esto es mejor que la multiplicación incluso para matrices medianas grandes, pero para un 2x2 o un 3x3, la sobrecarga adicional es un asesino. (No voy a adivinar dónde cae el punto de equilibrio, pero uno podría probarlo con la suficiente facilidad.)

+0

Prefiero la opción 1 o incluso la opción de bucle 2. Me sorprendería si la pérdida de velocidad del intérprete en el ciclo es significativa y prefiero el código más flexible. Imagine hacer un determinante de 5x5 de esta manera –

+2

@ChrisA - ¿Es 12.041883/0.002141 = 5624.4 una diferencia impactante? –

+0

+1 ¡Sí! Buena respuesta. Tal vez hay una forma más flexible de hacer esto ... ya que se pone realmente desagradable ya que la cantidad de términos aumenta con 'p!' –

3

me gustaría utilizar arrayfun:

d = arrayfun(@(w) det(Y(:, :, w)), 1 : size(Y, 3)); 

Editar: prueba de velocidad:

p = 10; 
n = 1e4; 
Y = rand(p,p,n); 

Prueba 1:

>> tic, d1 = arrayfun(@(w) det(Y(:, :, w)), 1 : size(Y, 3)); toc 
Elapsed time is 0.139030 seconds. 

Prueba 2 (por astillas de madera):

>> tic, X = squeeze(mat2cell(Y,p,p,ones(1,n))); d2= cellfun(@det,X); toc 
Elapsed time is 1.318396 seconds. 

Prueba 3 (enfoque ingenuo):

>> p = 10; 
>> n = 1e4; 
>> Y = rand(p,p,n); 
>> tic; d = nan(n, 1); for w = 1 : length(d), d(w) = det(Y(:, :, w)); end; toc 
Elapsed time is 0.069279 seconds. 

Conclusión: El enfoque ingenuo es el más rápido.

+0

¿Se puede usar el caso de prueba en @woodchips answer para hacer una comparación de velocidad entre esta sintaxis y las opciones cellfun y special-case? – tmpearce

+0

@tmpearce, no hay problema, agregué el resultado de la prueba – Serg

+0

Los resultados dependerán del tamaño de 'p': con p = 2, el modo @woodchips es probablemente bastante rápido. Sin embargo, es interesante que con p = 10, el ciclo es más rápido que arrayfun. – tmpearce