Esta va a ser una entrada «utilitaria»: poca explicación, pero aportando código fuente, por si alguien lo necesita. ¿Recuerda la entrada en la que presenté la descomposición de Cholesky? En aquel momento, no incluí el algoritmo de la descomposición, porque quería experimentar un poco con la implementación. Ya lo he hecho, y estoy más o menos satisfecho con el resultado.
Antes de ver el código, le explico el contexto. Este es un método de una estructura Matrix, que encapsula solamente una matriz bidimensional de valores flotantes (campo values). Dentro de la estructura se definen todos los operadores y métodos que podéis imaginar. En paralelo, defino otra estructura, LowerMatrix, para representar matrices triangulares inferiores. No hay relaciones de herencia, al tratarse de estructura, pero la clase de matrices triangulares permite definir métodos y operadores más eficientes. Es importante tener en cuenta que LowerMatrix gasta exactamente la misma memoria que Matrix. La ventaja está en esos métodos que evitan procesar las dos mitades de la matriz.
También hace falta saber que he definido operadores de conversión implícitos que transforman un array bidimensional en uno u otro tipo de matrices. Este operador es el que me permite evitar un constructor explícito para devolver un valor en el método:
public unsafe LowerMatrix Cholesky()
{
int n = Rows;
double[,] dest = new double[n, n];
double[,] src = values;
double* tmp = stackalloc double[n + n];
// First column is special.
double ajj = src[0, 0];
if (ajj <= 0)
{
dest[0, 0] = double.NaN;
return dest;
}
dest[0, 0] = ajj = Math.Sqrt(ajj);
double r = 1 / ajj;
n--;
for (int i = 1; i <= n; i++)
dest[i, 0] = src[i, 0] * r;
for (int j = 1; j <= n; j++)
{
// Compute the diagonal cell.
double v = 0.0;
for (int i = 0; i < j; i++)
{
double a = dest[j, i];
v += a * a;
}
ajj = src[j, j] - v;
if (ajj <= 0)
{
dest[j, j] = double.NaN;
return dest;
}
dest[j, j] = ajj = Math.Sqrt(ajj);
// Compute the other cells of column J.
if (j < n)
{
r = 1 / ajj;
for (int i = 0; i < j; i++)
tmp[i] = dest[j, i];
for (int i = j; i < n; i++)
{
v = 0.0;
for (int k = 0; k < j; k++)
v += dest[i + 1, k] * tmp[k];
tmp[i] = v;
}
for (int i = j; i < n; i++)
dest[i + 1, j] = (src[i + 1, j] - tmp[i]) * r;
}
}
return dest;
}
Sobre el algoritmo, en sí: el algoritmo de Cholesky puede fallar cuando la matriz de origen no es positiva semidefinida. Mi implementación detecta ese caso al calcular las raíces cuadradas… y se limita a parar, poniendo un NaN en la celda diagonal donde se ha detectado el problema. Esto significa que el método, en la práctica, asume que la matriz es positiva semidefinida. En mi código, tengo un segundo método, TryCholesky, que devuelve un valor lógico para ver si la conversión fue posible, y retorna la matriz transformada como parámetro de salida.
Desde el punto de vista de un programador de C#, el único detalle interesante es el uso de stackalloc para reservar un array de memoria en la pila, en vez de usar memoria dinámica. Esto es lo que obliga a declarar el método con unsafe.
En rendimiento, el método es más rápido que la «versión base» de la librería que he visto que es más rápida usando sólo C# (usando Intel MKL, no hay competencia posible). Me refiero a la versión base porque, para matrices grandes, las librerías serias suelen dividir la matriz en bloques que se pueden procesar en paralelo. Este código, como ve, no usa threads, instrucciones SIMD y sólo utiliza punteros para la caché. En menos palabras: todo es mejorable.