Tengo vectores de fila (densos y grandes) v_1, ..., v_D, y necesito que una matriz X sea diagonal de bloque (cuyo i-ésimo bloque es v_i), para realizar productos matriciales como

X.transpose() * H.selfadjointView<Lower>() * X  

Es decir, X debería ser escasa y su i-ésima fila es v_i:

(pseudo code)
RowVectorXd v1(1,2,3), v2(4,5), v3(6,7,8,9);
SparseMatrix<double> X(3,9); 

// I need X to be
X = 1 2 3 . . . . . .
    . . . 4 5 . . . .
    . . . . . 6 7 8 9 

// where . means 0

EDITAR: Mi pregunta es: ¿es posible hacer este producto sin formar X ni copiar las v_i, sino simplemente usando referencias a ellas? Mi preocupación es el rendimiento y la huella de memoria.

Creo que la solución tiene algo que ver con Eigen :: Map <... Stride ...> pero no puedo conseguirlo.

Muchas gracias.

1
itQ 16 ene. 2017 a las 22:28
Su pregunta sería más fácil de responder para la gente si incluyera un "ejemplo mínimo" del tipo de matriz que desea formar y el código que probó que sí usa la copia. ¡Buena suerte!
 – 
kcrisman
16 ene. 2017 a las 22:35
Gracias por tu comentario. Agregué el pseudocódigo correspondiente.
 – 
itQ
16 ene. 2017 a las 22:55
No puede copiar los vectores en SparseMatrix sin copiarlos (obviamente ...). Si su SparseMatrix tiene una estructura de bloque denso puro (es decir, sin ceros estructurales dentro de los bloques), puede acceder a estos bloques usando un Eigen::Map sin copiarlos. Si realmente necesita construir su matriz X depende de lo que pretenda hacer con ella.
 – 
chtz
17 ene. 2017 a las 12:10
Ese es exactamente el punto. Me pregunto si debería formarse X. Edité mi pregunta en consecuencia. Gracias.
 – 
itQ
17 ene. 2017 a las 14:01
Supongo que la pregunta sigue siendo, ¿qué pretendes hacer con X^T*H*X? Consistirá en subbloques de rango 1 (y almacenar explícitamente matrices de rango 1 a menudo es una mala idea, especialmente si conocía los vectores para construirlo). Además, ¿H es escaso o denso?
 – 
chtz
18 ene. 2017 a las 12:22

1 respuesta

La mejor respuesta

Terminé encontrando una respuesta a mi pregunta. Lo publico aquí en caso de que alguien más lo necesite. Para esta configuración particular de la matriz X, podemos formar directamente

B = X^T * H * X

Donde X es como en la pregunta y H es simétrico, sin formar explícitamente X. Como X ^ T H X es simétrico, solo necesitamos formar la diagonal inferior y usar

B.selfadjointView<Lower>() 

Llamar a B.

Tenga en cuenta que X ^ T H X es equivalente a H ** W, donde W es la siguiente matriz triangular inferior con bloques

W = v_1^T * v_1
    v_2^T * v_1    v_2^T * v_2
    .              .             v_3^T * v_3
    .              .             .               ...
    .              .             .               ...
    v_D^T * v_1    v_D^T * v_2   v_D^T * v_3     ...  v_D^T * v_D

Y ** es un operador sobrecargado que calcula el producto por componentes entre el escalar (i, j) de H y el bloque (i, j) en W. Entonces

B = H ** (v^T * v) 
0
itQ 23 ene. 2017 a las 13:37