Categorías
C#

Austra PolySolve

C’est la vie: elle est dure et souvent courte.

Es imposible escribir sobre ecuaciones algebraicas sin mencionar a Évariste Galois, quien no sólo cerró una larga historia de intentos de solución de este tipo de ecuaciones, sino que además tuvo una vida corta y trágica.

Dicen que los egipcios y los babilonios eran capaces de resolver las ecuaciones de segundo grado y, por supuesto, las lineales. Las de tercer y cuarto grado tuvieron que esperar al Renacimiento Italiano. Y luego, la teoría se estancó: nadie era capaz de resolver una ecuación de quinto grado general; sólo algunas versiones restringidas.

Lagrange estuvo a punto de demostrar que las ecuaciones de quinto grado y superiores no tenían una solución general. Fue Ruffini quien lo consiguió, aunque con algunos pequeños errores, que más tarde corrigió Abel. De todos modos, la teoría que Galois puso por escrito en 1830, cuando sólo tenía 18 años, tenía un alcance mucho mayor, y ofrecía una estructura más completa y versátil para estudiar las ecuaciones algebraicas. Mientras que el teorema de Ruffini-Abel se centraba en la solubilidad de ecuaciones algebraicas por medio de funciones elementales (como exponentes, logaritmos, etc.), la criatura que inventó nuestro héroe introdujo los llamados grupos de Galois para analizar la solubilidad y las simetrías en un marco más general. No solo abordó la solubilidad de ecuaciones, sino que la teoría es aplicable también en otras áreas de las matemáticas, como la teoría de números y la geometría algebraica.

Por desgracia, el artículo presentado en 1830 por nuestro héroe no tuvo el éxito que merecía. Cauchy le pasó la patata caliente a Poisson, y Poisson no entendió ni papa del tema. El rechazo cabreó a Galois, pero aprovechó para rescribir la demostración, y fue esta modificación la que finalmente fue reconocida, entre otros, por el propio Cauchy. Eso sí: tras la muerte de Galois…

Galois tenía una cabecita muy loca, le pirraba la política, y para colmo, estaba algo deprimido por la muerte de su padre. 1832 fue un año difícil para el chico. Estuvo en prisión un par de veces, por ciscarse en Louis Philippe, el penúltimo rey de Francia. Al salir de la cárcel por segunda vez, se enredó en un duelo absurdo, teóricamente por una coquette, aunque no es descartable que todo fuese una trampa de sus enemigos políticos. Se pasó la noche anterior al duelo escribiendo una carta sobre sus últimos avances matemáticos. Al día siguiente, interpuso su abdomen en la trayectoria de una bala, y sus contrincantes lo dejaron tirado sobre la hierba como a un chien. Un transeúnte lo vio y lo llevó al hospital, pero al día siguiente se reunió con su Creador, probablemente por culpa de una peritonitis.

And all the king’s horses and all the king’s men, couldn’t put Évariste together again.

Raíces reales

Mi curiosidad por estos temas viene de cuando tenía unos diez u once años: encontré la solución razonada de las ecuaciones de segundo grado, en un libro de electrónica, y me dio por intentar resolver por mi cuenta el problema de las cúbicas. No lo conseguí. Tropecé por casualidad con la sustitución de Vieta, pero no conseguí algo mucho más sencillo: cómo eliminar el término cuadrático, que suele ser el primer paso de la solución. Pero compré un libro que explicaba la fórmula cúbica y la cuártica, y me convertí en un friqui de las mates.

Volví a enredar con ecuaciones algebraicas en 2005, cuando me dio por probar si se podía escribir un ray tracer decente en C#. Es bastante frecuente tener que resolver ecuaciones de tercer y cuarto grado para calcular intersecciones entre rayos de luz y determinados tipos de objetos. La particularidad es que, en este contexto, sólo se necesitan las soluciones reales. Cuando las cosas se ponen feas, existe una técnica para encontrar las raíces reales de cualquier polinomio utilizando las secuencias de Sturm. Naturalmente, este algoritmo es una sólo una aproximación iterativa.

Raíces complejas, todas

Cuando estás escribiendo una librería como AUSTRA, te interesa resolver el problema más general, que es encontrar todas las raíces, ya sean complejas o reales, de un polinomio arbitrario. ¿Se acuerda de los valores propios? El método que utilizo en AUSTRA está basado en ellos.

Supongamos que queremos resolver la ecuación:

$$c_0 + c_1x + c_2x^2 + \cdots + c_{n-1}x^{n-1} + x^n$$El término de mayor grado está normalizado para que su coeficiente sea la unidad. Ahora formamos la siguiente matriz, conocida como «matriz de Frobenius»:

$$F=\pmatrix{0&0&0&\cdots&0&-c_0\cr
1&0&0&\cdots&0&-c_1\cr
0&1&0&\cdots&0&-c_2\cr
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\cr
0&0&0&\cdots&1&-c_{n-1}}$$Nos planteamos entonces encontrar los valores propios de $F$, que deben cumplir esta igualdad:

$$F\vec{v} = \lambda\vec{v}$$donde $\vec{v}$ es uno de los vectores propios. Si reordenamos los términos, nos encontramos con esto:

$$(F-\lambda I)\vec{v}=0$$donde $I$ es la matriz identidad. Para que esta igualdad se cumpla, el determinante de $(F-\lambda I)$ debe ser igual a cero. Y resulta que el determinante de $(F-\lambda I)$ es, precisamente, la ecuación original. Qué listo era Frobenius.

AUSTRA tiene un método muy eficiente para calcular valores propios, incluso en casos como estos, en los que la matriz no es simétrica. Por lo tanto, para resolver un polinomio primero lo normalizamos, luego creamos su matriz de Frobenius, y finalmente calculamos sus valores propios. La función global polySolve es la que se encarga de la implementación, en el lenguaje funcional de AUSTRA. En la aplicación de consola, podemos teclear lo siguiente:

> set v = [5, 4, 3, 2, 1]
ans ∊ ℝ(5)
5  4  3  2  1
> polysolve(v)
ans ∊ ℂ(4)
<0,137832; 0,678154>   <-0,537832; 0,358285>
<0,137832; -0,678154>  <-0,537832; -0,358285>

polySolve puede recibir tanto un vector con los coeficientes, como los coeficientes sueltos. En este caso, estamos resolviendo la ecuación de cuarto grado $5x^4+4x^3+3x^2+2x+1=0$, y el resultado son cuatro números complejos, conjugados a pares.

¿Quiere comprobar que las raíces son realmente soluciones de la ecuación? Hagamos esto entonces:

> polysolve(v).map(c => polyeval(c, v))
ans ∊ ℂ(4)
<-1,33227E-15; -7,77156E-16>   <-1,33227E-15; 4,44089E-16>
 <-1,33227E-15; 7,77156E-16>  <-1,33227E-15; -4,44089E-16>

polyEval sirve para evaluar un polinomio para un argumento complejo o real, y el método map crea un nuevo vector complejo calculando sus entradas con una función lambda, al estilo del método Select de LINQ. Incluso tenemos una función poliDerivative que, con los mismos argumentos que polyEval, evalúa la derivada del polinomio que le pasamos en la coordenada que le digamos. Esto, a su vez, es muy conveniente para buscar raíces reales con el método de Newton-Raphson… que también ofrece AUSTRA (función solve, a secas).

¿Librería o lenguaje?

Por supuesto, todo esto sería igual de fácil, eficiente y elegante, o quizás un poco más, si simplemente enchufásemos el package Austra.Library a un proyecto en .NET Core y utilizásemos directamente las clases. Pero he querido mostrar este ejemplo en el lenguaje de fórmulas de AUSTRA como demostración de un caso de uso importante para el lenguaje: es una forma rápida y sencilla de poner a prueba la funcionalidad de la librería.

Y hay más casos de uso, que explicaré más adelante.

Categorías
C#

El Gran Secreto de los Complejos

Al grano: el Gran Secreto de los Números Complejos es que, si quieres utilizar instrucciones AVX para acelerar los cálculos, la mejor forma de representarlos no es la que todos imaginamos: la parte real y, a continuación, la parte imaginaria.

Partamos de una regla básica de las instrucciones vectoriales:

  • Es mejor manejar estructuras de arrays que arrays de estructuras.

Observe que ésta es una píldora difícil de tragar en la Programación Orientada a Objetos. Complex es una clase que ya está (bien) definida en System.Numerics, pero para simplificar la explicación, voy a fingir que la definimos nosotros. Con la POO en mente, comenzaríamos definiendo la estructura, junto con sus métodos, y la haríamos probablemente implementar algunas interfaces, por completitud, en este plan:

public readonly struct Complex
{
    public double Real { get; }
    public double Imaginary { get; }

    public Complex(double re, double im) =>
        (Real, Imaginary) = (re, im);

    // Y así, sucesivamente…
}

Si quisiéramos entonces un vector de números complejos, haríamos algo parecido a esto:

public readonly struct ComplexVector
{
    private readonly Complex[] values;

    public unsafe ComplexVector(Complex[] values) =>
        this.values = values;

    // Y así, sucesivamente…
}

Pues bien, ahora llego yo (o la sacrosanta Realidad, si lo prefiere hacer menos personal) y le digo que la mejor forma de programar un vector de complejos, al menos si queremos acelerarlo con AVX, es la siguiente:

public readonly struct ComplexVector
{
    private readonly double[] re;
    private readonly double[] im;

    // Omito verificaciones de igual longitud
    // para simplificar el ejemplo.
    public ComplexVector(double[] re, double[] im) =>
        (this.re, this.im) = (re, im);

    public unsafe ComplexVector(Complex[] values)
    {
        this.re = new double[values.Length];
        this.im = new double[values.Length];
        fixed (double* p = re, q = im)
            for (int i = 0; i < values.Length; i++)
                (p[i], q[i]) = values[i];
    }

    // Y así, sucesivamente…
}

Podemos dejar la estructura Complex original: nos es útil. Pero al representar la lista de complejos, es mejor que cada campo vaya en su propia lista. Es cierto también que deberíamos utilizar AVX para convertir un array tradicional de complejos en un vector: existe la posibilidad, pero no lo voy a mostrar aquí, para simplificar. He omitido también un método de extensión que he añadido en una clase estática para poder "deconstruir" fácilmente un complejo en sus componentes. No tiene mucha trascendencia, pero ahí va, para que no haya tantos espacios en blanco:

[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static void Deconstruct(
    this Complex cmp, out double re, out double im) =>
    (re, im) = (cmp.Real, cmp.Imaginary);

¿Cómo se da cuenta uno?

¿Cómo se da cuenta uno de estas cosas? StackOverflow está lleno de consejos de este tipo, escritos por personas que ya lo han sufrido en sus carnes. Pero a uno no se le enciende la bombilla hasta que se pega con su propia pared. En este caso, fue intentando acelerar una Transformada Discreta de Fourier para AUSTRA. El código sin acelerar usaba los complejos como se han usado casi toda la vida: cada real con su parte imaginaria. A priori, cuando uno no conoce bien AVX, se imagina que será relativamente sencillo manejar dos números complejos dentro de un vector de cuatro valores reales, y que el conjunto de instrucciones va a estar de tu parte. Craso error.

Al final, tuve que limitarme a acelerar algunas partes que, oportunistamente, "se dejaban", casi siempre con código SSE y vectores de sólo 128 bits. Lo normal, cuando he acelerado otros algoritmos, ha sido reducir los tiempos de ejecuciones de cuatro hasta incluso ocho o diez veces. En este caso, la mejora sólo ha sido la mitad, en la mayoría de los casos, y en algunos, la tercera parte. Como resultado, tengo pendiente replantearme todo el asunto de la Transformada Discreta de Fourier, pero utilizando listas paralelas para los reales y los imaginarios.

Quiero que vea, de todas maneras, lo sencillo que es usar la técnica de estructura de listas en vez de listas de estructuras. El siguiente método es el producto escalar de dos vectores complejos, representados como Dios manda:

public static unsafe Complex operator *(
    ComplexVector v1, ComplexVector v2)
{
    if (v1.Length != v2.Length)
        throw new VectorLengthException();
    fixed (double* pr = v1.re, pi = v1.im,
                   qr = v2.re, qi = v2.im)
    {
        double sumRe = 0, sumIm = 0;
        int i = 0, size = v1.Length;
        if (Avx.IsSupported)
        {
            Vector256<double> accRe = Vector256<double>.Zero;
            Vector256<double> accIm = Vector256<double>.Zero;
            for (int top = size & ~3; i < top; i += 4)
            {
                var vpr = Avx.LoadVector256(pr + i);
                var vpi = Avx.LoadVector256(pi + i);
                var vqr = Avx.LoadVector256(qr + i);
                var vqi = Avx.LoadVector256(qi + i);
                accRe = Avx.Add(accRe,
                    Avx.Multiply(vpr, vqr)
                       .MultiplyAdd(vpi, vqi));
                accIm = Avx.Add(accIm,
                    Avx.Multiply(vpi, vqr)
                       .MultiplySub(vpr, vqi));
            }
            sumRe = accRe.Sum();
            sumIm = accIm.Sum();
        }
        for (; i < size; i++)
        {
            sumRe += pr[i] * qr[i] + pi[i] * qi[i];
            sumIm += pi[i] * qr[i] - pr[i] * qi[i];
        }
        return new(sumRe, sumIm);
    }
}

Si no lo recuerda del álgebra, cuando se trata de vectores complejos, el producto escalar usa la conjugada del segundo operando. Esto es: la parte imaginaria del segundo operando invierte su signo. Esto es lo que permite que el producto escalar de un vector consigo mismo sea un valor real.

El código vectorial tiene una correspondencia uno a uno con el código escalar que maneja la parte de los arrays que puede sobrar al final. Si repasa la entrada anterior, la del algoritmo de Welford, notará el mismo patrón: a pesar de que el algoritmo escalar es bastante oscuro, es relativamente sencillo convertir esa parte en código vectorial. La parte más complicada de la entrada pasada era cuando había que mezclar los cuatro acumuladores individuales. Es el mismo problema que hemos visto ya unas cuantas veces cuando calculamos un producto escalar: acumular es sencillo. Lo complicado es sumar después los cuatro acumuladores.

El Misterioso Constructor Postergado

Para no dejar demasiados cabos sueltos, aquí tiene una posible implementación del código que reparte partes reales e imaginarias a sus respectivos arrays. Mi primer impulso fue utilizar Avx2.GatherVector: leer cuatro partes reales saltándome las imaginarias, y luego leer cuatro partes imaginarias. Pero, por desgracia, el tiempo de ejecución del constructor se disparó al doble. No hay forma humana de predecir estas cosas, que no sea prueba, benchmark y error.

La versión que funciona, y que reduce casi a la mitad el tiempo de la versión de más arriba, lee cuatro complejos en dos vectores de 256 bits. Lo primero que hace es usar Avx.Shuffle para "barajar las cartas" y juntar todas las partes reales en un mismo vector de 256 bits, y las imaginarias en otro. No me pida que calcule estas cosas de memorias. Cuando me tocan estos marrones, tengo todavía que pillar un cuaderno y lápiz, e irme a páginas como ésta, para repasar los diagramas. He visto también que el Shuffle se puede sustituir también por llamadas a UnpackHigh/UnpackLoad. Es probable que estas llamadas den mejores tiempos, pero no me ha dado tiempo a hacer la prueba.

El problema de Shuffle (y de las alternativas mencionadas) es que te deja los números en el orden [1, 3, 2, 4]. Si no es importante respetar el orden, se pueden quedar así. Pero si hay que reordenar los elementos, hay que usar Avx2.Permute4x64 para ello. En general, AVX intenta, dentro de lo posible, no pasar valores de una mitad a la otra mitad del vector. Hay que usar cosas introducidas en AVX2 para conseguirlo. Por ese motivo, el constructor verifica si Avx2.IsSupported antes de lanzarse al río:

public unsafe ComplexVector(Complex[] values)
    : this(values.Length)
{
    fixed (double* p = re, q = im)
    fixed (Complex* r = values)
    {
        int i = 0;
        if (Avx2.IsSupported)
        {
            for (int top = values.Length & ~7; i < top; i += 4)
            {
                var v1 = Avx.LoadVector256((double*)(r+i));
                var v2 = Avx.LoadVector256((double*)(r+i+2));
                Avx.Store(p + i, Avx2.Permute4x64(
                    Avx.Shuffle(v1, v2, 0b0000), 0b11011000));
                Avx.Store(q + i, Avx2.Permute4x64(
                    Avx.Shuffle(v1, v2, 0b1111), 0b11011000));
            }
        }
        for (; i < values.Length; i++)
            (p[i], q[i]) = r[i];
    }
}

Números: en mi máquina, un i9-11900K, crear un ComplexVector directamente a partir de un array de 1024 complejos, tardaba más o menos un milisegundo. Con las mejoras AVX2, tarda 650 microsegundos. Casi la mitad. Y lo mejor, para mi gusto, es que no he tenido que usar paralelismo con tareas. El usuario de la librería ya usará ese paralelismo cuando lo considere necesario, y tendrá las manos más libres.

Como regalo, le dejo la conversión inversa: de vector complejo a array de complejos:

public unsafe static explicit
    operator Complex[](ComplexVector v)
{
    Complex[] result = new Complex[v.Length];
    fixed (double* p = v.re, q = v.im)
    fixed (Complex* r = result)
    {
        int i = 0;
        if (Avx2.IsSupported)
        {
            for (int top = v.Length & ~3; i < top; i += 4)
            {
                var vr = Avx.LoadVector256(p + i);
                var vi = Avx.LoadVector256(q + i);
                Avx.Store((double*)(r + i),
                    Avx2.Permute4x64(Avx.Permute2x128(
                    vr, vi, 0b0010_0000), 0b11_01_10_00));
                Avx.Store((double*)(r + i + 2),
                    Avx2.Permute4x64(Avx.Permute2x128(
                    vr, vi, 0b0011_0001), 0b11_01_10_00));
                }
            }
        for (; i < result.Length; i++)
            r[i] = new(p[i], q[i]);
    }
    return result;
}

El tiempo de ejecución se reduce "sólo" a las tres cuartas partes, pero yo creo que merece la pena.