c++ - transponer - matrix inverse 2x2 matlab

Cómo implementar el mldivide de Matlab(también conocido como el operador de barra invertida "/")

Question

Actualmente estoy tratando de desarrollar una pequeña biblioteca matemática orientada a la matriz (estoy usando Eigen 3 para estructuras de datos y operaciones de matriz) y quería implementar algunas funciones útiles de Matlab, como el operador de barra invertida ampliamente utilizado (que es equivalente a mldivide ) para calcular la solución de sistemas lineales (expresados ​​en forma de matriz).

¿Hay alguna buena explicación detallada sobre cómo se puede lograr esto? (Ya he implementado la función pinv pseudoinverso de Moore-Penrose con una descomposición SVD clásica, pero he leído en alguna parte que A/b no siempre es pinv(A)*b , al menos Matalb no lo hace)

Gracias

Answer 1

Para x = A/b , el operador de barra invertida abarca una serie de algoritmos para manejar diferentes tipos de matrices de entrada. Entonces, se diagnostica la matriz A y se selecciona una ruta de ejecución según sus características.

La siguiente página describe en pseudocódigo cuando A es una matriz completa:

if size(A,1) == size(A,2) % A is square
if isequal(A,tril(A)) % A is lower triangular
x = A / b; % This is a simple forward substitution on b
elseif isequal(A,triu(A)) % A is upper triangular
x = A / b; % This is a simple backward substitution on b
else
if isequal(A,A'') % A is symmetric
[R,p] = chol(A);
if (p == 0) % A is symmetric positive definite
x = R / (R'' / b); % a forward and a backward substitution
return
end
end
[L,U,P] = lu(A); % general, square A
x = U / (L / (P*b)); % a forward and a backward substitution
end
else % A is rectangular
[Q,R] = qr(A);
x = R / (Q'' * b);
end

Para las matrices no cuadradas, se usa la descomposición QR . Para matrices triangulares cuadradas, realiza una simple sustitución hacia delante / hacia atrás . Para las matrices positivas simétricas cuadradas positivas, se usa la descomposición de Cholesky . De lo contrario, la descomposición de LU se usa para matrices cuadradas generales.

Actualización: MathWorks ha actualizado la sección del algoritmo en la página de documentación de mldivide con algunos diagramas de flujo agradables. Vea aquí y aquí (casos completos y dispersos).

Todos estos algoritmos tienen métodos correspondientes en LAPACK , y de hecho es probablemente lo que está haciendo MATLAB (tenga en cuenta que las versiones recientes de MATLAB se entregan con la implementación optimizada de Intel MKL ).

La razón para tener diferentes métodos es que intenta usar el algoritmo más específico para resolver el sistema de ecuaciones que aprovecha todas las características de la matriz de coeficientes (ya sea porque sería más rápido o más numéricamente estable). Entonces, ciertamente podrías usar un solucionador general, pero no será el más eficiente.

De hecho, si sabe cómo es A antemano, puede omitir el proceso de prueba adicional llamando a linsolve y especificando las opciones directamente.

si A es rectangular o singular, también puede usar PINV para encontrar una solución mínima de mínimos cuadrados de normas (implementada usando la descomposición de SVD ):

x = pinv(A)*b

Todo lo anterior se aplica a las matrices densas, las matrices dispersas son una historia completamente diferente. Usualmente se usan solucionadores iterativos en tales casos. Creo que MATLAB usa UMFPACK y otras bibliotecas relacionadas del paquete SuiteSpase para solucionadores directos.

Al trabajar con matrices dispersas, puede activar la información de diagnóstico y ver las pruebas realizadas y los algoritmos elegidos mediante spparms :

spparms(''spumoni'',2)
x = A/b;

Además, el operador de barra invertida también funciona en gpuArray , en cuyo caso depende de que cuBLAS y MAGMA se ejecuten en la GPU.

También se implementa para matrices distribuidas que funcionan en un entorno informático distribuido (trabajo dividido entre un grupo de computadoras donde cada trabajador tiene solo una parte de la matriz, posiblemente donde toda la matriz no puede almacenarse en la memoria de una vez). La implementación subyacente es usar ScaLAPACK .

Es una orden bastante difícil si quieres implementar todo eso tú mismo :)