Optimización Super C ++ de multiplicación de matrices con Armadillo

Estoy usando Armadillo para hacer multiplicaciones matriciales muy intensas con longitudes laterales , donde puede ser de hasta 20 o incluso más. Estoy usando Armadillo con OpenBLAS para la multiplicación de matrices, que parece estar haciendo un muy buen trabajo en núcleos paralelos, excepto que tengo un problema con el formalismo de la multiplicación en Armadillo para una súper optimización del rendimiento.$2^n$$n$

Digamos que tengo un bucle en la siguiente forma:

arma::cx_mat stateMatrix, evolutionMatrix; //armadillo complex matrix type
for(double t = t0; t < t1; t += 1/sampleRate)
{
    ...
    stateMatrix = evolutionMatrix*stateMatrix;
    ...
}

En C ++ fundamental, el problema aquí es que C ++ asignará un nuevo objeto de cx_matalmacenar evolutionMatrix*stateMatrixy luego copiará el nuevo objeto stateMatrixcon operator=(). Esto es muy, muy ineficiente. Es bien sabido que devolver clases complejas de tipos de datos grandes es una mala idea, ¿verdad?

La forma en que veo que esto es más eficiente es con una función que realiza la multiplicación en la forma:

void multiply(const cx_mat& mat1, const cx_mat& mat2, cx_mat& output)
{
    ... //multiplication of mat1 and mat2 and then store it in output
}

De esta manera, uno no tiene que copiar objetos enormes con valor de retorno, y la salida no tiene que reasignarse con cada multiplicación.

La pregunta : ¿Cómo puedo encontrar un compromiso, en el que pueda usar Armadillo para la multiplicación con su agradable interfaz de BLAS, y hacer esto de manera eficiente sin tener que recrear objetos de matriz y copiarlos con cada operación?

¿No es este un problema de implementación en Armadillo?

En C ++ fundamental, el problema aquí es que C ++ asignará un nuevo objeto de cx_mat para almacenar evolutionMatrix * stateMatrix, y luego copiará el nuevo objeto a stateMatrix con operator = ().

Creo que tienes razón en que está creando temporarios, que es demasiado lento, pero creo que la razón por la que está haciendo eso es incorrecta.

Armadillo, como cualquier buena biblioteca de álgebra lineal de C ++, usa plantillas de expresión para implementar la evaluación retardada de expresiones. Al escribir un producto de matriz como A*B, no se crean los temporales, en lugar del armadillo hace que un objeto temporal ( x) que mantiene las referencias a Ay B, a continuación, más tarde, dada una expresión como C = x, calcula el producto de la matriz que almacena el resultado directamente en C, sin crear ningún temporales.

También utiliza esta optimización para manejar expresiones como A*B*C*D, donde, dependiendo del tamaño de la matriz, ciertos órdenes de multiplicación son más eficientes que otros.

¿No es este un problema de implementación en Armadillo?

Si Armadillo no está realizando esta optimización, sería un error en Armadillo que debería informarse a los desarrolladores.

Sin embargo, en su caso, hay otro problema que es más importante. En una expresión como A=B*Cel almacenamiento de Ano contiene ningún dato de entrada si Ano tiene alias Bo C. En su caso, A = A*Bescribir cualquier cosa en la matriz de salida también modificaría una de las matrices de entrada.

Incluso dada su función sugerida

multiply(const cx_mat&, const cx_mat&, cx_mat&)

¿Cómo exactamente ayudaría esa función en la expresión multiply(A, B, A)? Para la mayoría de las implementaciones ordinarias de esa función, eso conduciría a un error. Tendría que usar algo de almacenamiento temporal por sí solo, para asegurarse de que sus datos de entrada no estén dañados. Su sugerencia es más o menos cómo Armadillo ya implementa la multiplicación de matrices, pero creo que probablemente tenga cuidado de evitar situaciones como la multiply(A, B, A)asignación de un temporal.

La explicación más probable de por qué Armadillo no está haciendo esta optimización es que sería incorrecto hacerlo.

Finalmente, hay una manera mucho más simple de hacer lo que quieres, como esta:

cx_mat *A, *Atemp, B;
for (;;) {
  *Atemp = (*A)*B;
  swap(A, Atemp);
}

Esto es idéntico a

cx_mat A, B;
for (;;) {
  A = A*B;
}

pero asigna una matriz temporal, en lugar de una matriz temporal por iteración.

@BillGreene apunta a la "optimización del valor de retorno" como una forma de evitar el problema fundamental, pero esto en realidad solo ayuda a la mitad. Suponga que tiene un código de esta forma:

struct ExpensiveObject { ExpensiveObject(); ~ExpensiveObject(); };

ExpensiveObject operator+ (ExpensiveObject &obj1,
                           ExpensiveObject &obj2)
{
   ExpensiveObject tmp;
   ...compute tmp based on obj1 and obj2...
   return tmp;
}

void f() {
  ExpensiveObject o1, o2, o3;
  ...initialize o1, o2...;
  o3 = o1 + o2;
}

Un ingenuo compilador

1. crear una ranura para almacenar el resultado de la operación plus (una temporal),
2. operador de llamada +,
3. cree el objeto 'tmp' dentro del operador + (un segundo temporal),
4. calcular tmp,
5. copia tmp en la ranura de resultados,
6. destruir tmp,
7. copie el objeto resultante en o3
8. destruir el objeto resultante

La optimización del valor de retorno solo puede unificar el objeto 'tmp' y la ranura 'resultado', pero no elimina la necesidad de una copia. Entonces, todavía te queda la creación de un temporal, la operación de copia y la destrucción de un temporal.

La única forma de evitar esto es que operator + no devuelve un objeto, sino un objeto de alguna clase intermedia que, cuando se asigna a un ExpensiveObject, realiza la operación de suma y copia en su lugar. Este es el enfoque típico utilizado en las bibliotecas de plantillas de expresión .

Stackoverflow ( https://stackoverflow.com/ ) es probablemente un mejor foro de discusión para esta pregunta. Sin embargo, aquí hay una respuesta corta.

Dudo que el compilador de C ++ esté generando código para esta expresión como usted describe anteriormente. Todos los compiladores modernos de C ++ implementan una optimización llamada "optimización del valor de retorno" ( http://en.wikipedia.org/wiki/Return_value_optimization ). Con la optimización del valor de retorno, el resultado de evolutionMatrix*stateMatrixse almacena directamente en stateMatrix; No se realiza ninguna copia.

Obviamente, existe una considerable confusión sobre este tema y esa es una de las razones por las que sugerí que Stackoverflow podría ser un mejor foro. Hay muchos "abogados de lenguaje" en C ++, mientras que la mayoría de nosotros preferimos pasar nuestro tiempo en CSE. ;-)

Creé el siguiente ejemplo simple basado en la publicación del profesor Bangerth:

#ifndef NDEBUG
#include <iostream>

using namespace std;
#endif

class ExpensiveObject  {
public:
  ExpensiveObject () {
#ifndef NDEBUG
    cout << "ExpensiveObject  constructor called." << endl;
#endif
    v = 0;
  }
  ExpensiveObject (int i) { 
#ifndef NDEBUG
    cout << "ExpensiveObject  constructor(int) called." << endl;
#endif
    v = i; 
  }
  ExpensiveObject (const ExpensiveObject  &a) {
    v = a.v;
#ifndef NDEBUG
    cout << "ExpensiveObject  copy constructor called." << endl;
#endif
  }
  ~ExpensiveObject() {
#ifndef NDEBUG
    cout << "ExpensiveObject  destructor called." << endl;
#endif
  }
  ExpensiveObject  operator=(const ExpensiveObject  &a) {
#ifndef NDEBUG
    cout << "ExpensiveObject  assignment operator called." << endl;
#endif
    if (this != &a) {
      return ExpensiveObject (a);
    }
  }
  void print() const {
#ifndef NDEBUG
    cout << "v=" << v << endl;
#endif
  }
  int getV() const {
    return v;
  }
private:
  int v;
};

ExpensiveObject  operator+(const ExpensiveObject  &a1, const ExpensiveObject  &a2) {
#ifndef NDEBUG
  cout << "ExpensiveObject  operator+ called." << endl;
#endif
  return ExpensiveObject (a1.getV() + a2.getV());
}

int main()
{
  ExpensiveObject  a(2), b(3);
  ExpensiveObject  c = a + b;
#ifndef NDEBUG
  c.print();
#endif
}

Parece más complicado de lo que realmente es porque quería eliminar completamente todo el código para imprimir la salida al compilar en modo optimizado. Cuando ejecuto la versión compilada con una opción de depuración, obtengo el siguiente resultado:

ExpensiveObject  constructor(int) called.
ExpensiveObject  constructor(int) called.
ExpensiveObject  operator+ called.
ExpensiveObject  constructor(int) called.
v=5
ExpensiveObject  destructor called.
ExpensiveObject  destructor called.
ExpensiveObject  destructor called.

Lo primero que debe notar es que no se construyen temporarios, solo a, by c. El constructor predeterminado y el operador de asignación nunca se llaman porque no son necesarios en este ejemplo.

El profesor Bangerth mencionó plantillas de expresión. De hecho, esta técnica de optimización es muy importante para obtener un buen rendimiento en una biblioteca de clases de matriz. Pero es importante solo cuando las expresiones de objeto son más complicadas que simplemente a + b. Si, por ejemplo, mi prueba fue en su lugar:

  ExpensiveObject  a(2), b(3), c(9);
  ExpensiveObject  d = a + b + c;

Obtendría el siguiente resultado:

ExpensiveObject  constructor(int) called.
 ExpensiveObject  constructor(int) called.
 ExpensiveObject  constructor(int) called.
 ExpensiveObject  operator+ called.
 ExpensiveObject  constructor(int) called.
 ExpensiveObject  operator+ called.
 ExpensiveObject  constructor(int) called.
 ExpensiveObject  destructor called.
 v=14
 ExpensiveObject  destructor called.
 ExpensiveObject  destructor called.
 ExpensiveObject  destructor called.
 ExpensiveObject  destructor called.

Este caso muestra la construcción indeseable de un temporal (5 llamadas al constructor y dos llamadas del operador +). El uso adecuado de las plantillas de expresión (un tema que está más allá del alcance de este foro) evitaría esto temporalmente. (Para los altamente motivados, se puede encontrar una discusión particularmente legible de las plantillas de expresión en el capítulo 18 de http://www.amazon.com/C-Templates-The-Complete-Guide/dp/0201734842 ).

Finalmente, la verdadera "prueba" de lo que realmente está haciendo el compilador proviene del examen del código de ensamblaje que genera el compilador. Para el primer ejemplo, cuando se compila en modo optimizado, este código es asombrosamente simple. Todas las llamadas a funciones se han optimizado y el código de ensamblado esencialmente carga 2 en un registro, 3 en un segundo y las agrega.

$\mathcal{O}(n^{2.8})$$\mathcal{O}(n^2)$$n$

Es decir, a menos que incurra en una enorme constante en la copia, que en realidad no es tan descabellada, porque la versión con copia es mucho más costosa en otro aspecto: necesita mucha más memoria. Entonces, si termina teniendo que cambiar hacia y desde el disco duro, la copia podría convertirse en el cuello de botella. Sin embargo, incluso si no copia nada usted mismo, un algoritmo fuertemente paralelizado puede hacer algunas copias por sí mismo. Realmente, la única forma de asegurarse de que no se utilizará demasiada memoria en cada paso es dividir la multiplicación en columnasstateMatrix , de modo que solo se hagan pequeñas multiplicaciones a la vez. Por ejemplo, puedes definir

class HChunkMatrix // optimised for destructive left-multiplication updates
{
  std::vector<arma::cx_mat> colChunks; // e.g. for an m×n matrix,
                                      //  use √n chunks, each an m×√n matrix
 public:
  ...

  HChunkMatrix& operator *= (const arma::cx_mat& lhMult) {
    for (&colChunk: colChunks) {
      colChunk = lhMult * colChunk;
    }
    return *this;
  }
}

También debe considerar si incluso necesita evolucionar eso stateMatrixcomo uno en primer lugar. Si básicamente solo desea una evolución temporal independiente de nlos kets de estado, entonces también podría evolucionarlos uno por uno, lo que es mucho menos costoso para la memoria. En particular, si evolutionMatrixes escaso , lo que definitivamente deberías revisar. Porque esto es básicamente un hamiltoniano, ¿no? Los hamiltonianos a menudo son escasos o aproximadamente escasos.

_{$\mathcal{O}(n^{2.38})$}

Modern C ++ tiene una solución para el problema mediante el uso de "mover constructores" y "referencias de valor".

Un "constructor de movimiento" es un constructor para una clase, por ejemplo, una clase de matriz, que toma otra instancia de la misma clase y mueve los datos de la otra instancia a la nueva instancia, dejando vacía la instancia original. Típicamente, un objeto matriz tendrá dos números para el tamaño y un puntero a los datos. Cuando un constructor normal duplicaría los datos, un constructor de movimiento solo copiará los dos números y el puntero, por lo que esto es realmente rápido.

Una "referencia de valor", escrita por ejemplo como "matriz &&" en lugar de la "matriz &" habitual se usa para variables temporales. Declararía una multiplicación matricial como la devolución de una matriz &&. Al hacerlo, el compilador se asegurará de que se utilizará un constructor de movimientos muy barato para obtener el resultado de la función que lo llama. Entonces, una expresión como resultado = (a + b) * (c + d) donde a, b, c, d son objetos de matriz enormes, sucederá sin ninguna copia.

Buscar en Google "referencias de valor y constructores de movimientos" encontrará ejemplos y tutoriales.

$v'*M*M*M*M*M$$M*M*M*M*M*v$

Por otra parte, creo que OpenBLAS tiene una colección más grande de optimizaciones específicas de arquitectura, por lo que Eigen puede o no ser una victoria para usted. Desafortunadamente, no hay una biblioteca de álgebra lineal tan impresionante que ni siquiera tenga que considerar a los demás cuando lucha por "el último 10%" de rendimiento. Los envoltorios no son una solución al 100%; la mayoría (¿todos?) de ellos no pueden aprovechar la capacidad de eigen para fusionar cálculos de esta manera.