c++ - transpuestas - matriz transpuesta

Tipo Eigen eficiente para matriz simétrica densa.

Question

¿Tiene Eigen un tipo eficiente para la matriz densa, de tamaño fijo, simétrica? (hey, son omnipresentes!)

Es decir, para N = 9, debe almacenar solo (1 + 9) * 9/2 == 45 elementos y tiene las operaciones apropiadas. Por ejemplo, debería haber una adición eficiente de dos matrices simétricas, que devuelve una matriz simétrica similar.

Si no hay tal cosa, ¿qué acciones (se parece a this ) debería hacer para introducir ese tipo en Eigen? ¿Tiene conceptos de "Vistas"? ¿Puedo escribir algo así como "vista de matriz" para mi propio tipo, que lo haría Eigen-friednly?

PS Probablemente puedo tratar la matriz plana como una matriz 1xN usando el map y hacer operaciones en ella. Pero no es la solución más limpia.

Answer 1

Sí, eigen3 tiene el concepto de views . Aunque no hace nada para el almacenamiento. Sin embargo, como idea, puede compartir un bloque más grande para dos matrices simétricas del mismo tipo:

Matrix<float,4,4> A1, A2; // assume A1 and A2 to be symmetric
Matrix<float,5,4> A;
A.topRightCorner<4,4>().triangularView<Upper>() = A1;
A.bottomLeftCorner<4,4>().triangularView<Lower>() = A2;

Aunque es bastante engorroso, y solo lo usaría si tu memoria es realmente preciosa.

Answer 2

El almacenamiento empaquetado de matrices simétricas es un gran enemigo del código vectorizado, es decir, de la velocidad. La práctica estándar es almacenar los coeficientes N * (N + 1) / 2 relevantes en la parte triangular superior o inferior de una matriz NxN de densidad completa y dejar el resto (N-1) * N / 2 sin referencia. Todas las operaciones en la matriz simétrica se definen teniendo en cuenta este almacenamiento peculiar. En eigen, usted tiene el concepto de vistas triangulares y auto-adjuntas para obtener esto.

A partir de la referencia eigen : (para matrices reales selfjojo == simétricas).

Al igual que para la matriz triangular, puede hacer referencia a cualquier parte triangular de una matriz cuadrada para verla como una matriz de autounión y realizar operaciones especiales y optimizadas. De nuevo, nunca se hace referencia a la parte triangular opuesta y se puede utilizar para almacenar otra información.

A menos que la memoria sea un gran problema, sugeriría dejar la parte sin referencia de la matriz vacía. (Código más legible, sin problemas de rendimiento).

Answer 3

Tipo eficiente para matriz simétrica.

Simplemente asigna valores a las partes triangulares inferior / superior de la matriz y utiliza las vistas Eigen triangular y autoadjoint. Sin embargo, he probado ambos en matrices pequeñas de tamaño fijo. Me di cuenta de que, en lo que respecta al rendimiento, utilizar vistas no siempre es la mejor opción Considere el siguiente código:

Eigen::Matrix2d m;
m(0,0) = 2.0;
m(1,0) = 1.0;
// m(0,1) = 1.0;
m(1,1) = 2.0;
Eigen::Vector2d v;
v << 1.0,1.0;
auto result = m.selfadjointView<Eigen::Lower>()*v;

El producto en la última línea es bastante lento en comparación con las soluciones alternativas que se presentan a continuación (aproximadamente un 20% más lento para matrices double 2x2 en mi caso). (El producto que usa la matriz completa, al descomentar m(0,1) = 1.0; y al usar auto result = m*v , es aún más rápido para matrices double 2x2 ).

EDIT : me olvidé de aliasing. auto result.noalise() acelera las cosas (las alternativas a continuación son aún más rápidas).

Algunas alternativas.

1) Almacenar matriz simétrica en vector.

Puede almacenar su matriz en un vector de tamaño 45. Sumar 2 matrices en formato vectorial es sencillo (solo sume los vectores). Pero tienes que escribir tu propia implementación para productos.

Aquí está la implementación de dicho producto de matrix * vector (denso, de tamaño fijo) donde la parte inferior de la matriz se almacena en forma de columna en un vector:

template <typename T, size_t S>
Eigen::Matrix<T,S,1> matrixVectorTimesVector(const Eigen::Matrix<T,S*(S+1)/2,1>& m, const Eigen::Matrix<T,S,1>& v)
{
    Eigen::Matrix<T,S,1> ret(Eigen::Matrix<T,S,1>::Zero());
    int counter(0);
    for (int i=0; i<S; ++i)
    {
        ret[i] += m(counter++)*v(i);
        for (int j=i+1; j<S; ++j)
        {
            ret[i] += m(counter)*v(j);
            ret[j] += m(counter++)*v(i);
        }
    }
    return ret;
}

2) Almacene solo la parte triangular e implemente sus propias operaciones.

Por supuesto, también podría implementar su propio matrix * vector producto, donde la matriz solo almacena los 45 elementos (supongamos que almacenamos la parte triangular inferior). Esta podría ser la solución más elegante, ya que mantiene el formato de una matriz (en lugar de usar un vector que representa una matriz). También puede usar funciones Eigen como en el siguiente ejemplo:

template <typename T, size_t S>
Eigen::Matrix<T,S,S> symmMatrixPlusSymmMatrix( Eigen::Matrix<T,S,S>& m1, const Eigen::Matrix<T,S,S>& m2)
{
    Eigen::Matrix<T,S,S> ret;
    ret.template triangularView<Eigen::Lower>() = m1 + m2; // no performance gap here!
    return ret;
}

En la función anterior (suma de 2 matrices simétricas), solo se visitan las partes triangulares inferiores de m1 y m2. Tenga en cuenta que la vista triangularView no presenta una brecha de rendimiento en este caso (afirmo esto en base a mis puntos de referencia).

En cuanto al producto matrix * vector , vea el ejemplo a continuación (el mismo rendimiento que el producto en la alternativa 1). El algoritmo solo visita la parte triangular inferior de la matriz.

template <typename T, size_t S>
Eigen::Matrix<T,S,1> symmMatrixTimesVector(const Eigen::Matrix<T,S,S>& m, const Eigen::Matrix<T,S,1>& v)
{
    Eigen::Matrix<T,S,1> ret(Eigen::Matrix<T,S,1>::Zero());
    int counter(0);

    for (int c=0; c<S; ++c)
    {
        ret(c) += m(c,c)*v(c);
        for (int r=c+1; r<S; ++r)
        {
            ret(c) += m(r,c)*v(r);
            ret(r) += m(r,c)*v(c);
        }
    }
    return ret;
}

La ganancia de rendimiento para el producto Matrix2d*Vector2d en comparación con el producto que utiliza la matriz completa (2x2 = 4 elementos) es del 10% en mi caso.