Julia

Multiple dispatch

Author

David Gómez-Castro

Diferenciación automática

Números duales

Al hacer un desarrollo de Taylor \[f(x+\delta) = f(x) + f'(x) \delta + \delta^2 \left( \frac{f''(x)}{2} + ... \right)\] Para entender la derivada querríamos poder “tirar” los términos superiores

Sea \(\varepsilon\) un elemento formal no-nulo tal que \(\varepsilon^2=0\).

Podemos construir los números duales \(a + b \varepsilon\) con \(a,b \in \mathbb R\) y \(\varepsilon \notin \mathbb R\).

Ejemplo tomado de Video de Youtube del Prof Alan Edelman (MIT)

Estructuras en Julia

Creamos una estructura de pares \((a,b)\) y la nombramos Dual

struct Dual <: Number
    f::Tuple{Float64,Float64}
end

El constructor se genera de manera automática

Dual( (1.0,2.1) )
Dual((1.0, 2.1))

Extendemos las funciones suma y resta a este tipo de elementos

import Base: + , -
+(x::Dual, y::Dual) = Dual( x.f .+ y.f );
-(x::Dual, y::Dual) = Dual( x.f .- y.f );

Y lo probamos

a = Dual( (1.0, 2.0) )
a+a
Dual((2.0, 4.0))

Bajo la premisa \(\varepsilon^2 = 0\), para multiplicar: \[(a+b\varepsilon)(c+d\varepsilon) = ac + (bc + ad)\varepsilon + bd \varepsilon^2 = ac + (bc + ad)\varepsilon.\]

Para dividir hay un truco sencillo \[ \frac{a+b\varepsilon}{c+d\varepsilon} = \frac{(a+b\varepsilon)(c-d\varepsilon)}{(c+d\varepsilon)(c-d\varepsilon)} = \frac{(a+b\varepsilon)(c-d\varepsilon)}{c^2-d^2 \varepsilon^2} = = \frac{(a+b\varepsilon)(c-d\varepsilon)}{c^2} \]

Finalmente \[ \frac{a+b\varepsilon}{c+d\varepsilon} = \frac{a}{c} + \frac{bc - ad}{c^2} \varepsilon \]

Nótese si \(f(x) = a + bx, g(x) = c + dx\) entonces \(\frac{d}{dx} \frac{f}{g} = \frac{bc - ad}{c^2}\).

import Base: * , /
*(x::Dual, y::Dual) = Dual( ( x.f[1]*y.f[1] , x.f[1]*y.f[2] + x.f[2]*y.f[1]) )   
/(x::Dual, y::Dual) = Dual(( 
        x.f[1] / y.f[1] , 
        (x.f[2]*y.f[1] - x.f[1]*y.f[2]) / (y.f[1]^2)
));

Pinceladas de Álgebra

Recordamos el conjunto de números complejos \(\mathbb C\) es

\(a + b i\) donde \(i \notin \mathbb R\)

y coincide con

\[ \mathbb R[x] / \langle x^2 + 1 \rangle \]

Del mismo modo, \(a + b \varepsilon\) donde \(\varepsilon \notin \mathbb R\) y \(\varepsilon^2 = 0\) es lo mismo que

\[ \mathbb R[x] / \langle x^2 \rangle \]

Un poco de maquillaje

Para verlo bonito

Base.show(io::IO,x::Dual) = print(io,x.f[1]," + ",x.f[2]," ε")
a
1.0 + 2.0 ε

Definimos una variable llamada \(\varepsilon\)

ϵ=Dual((0,1))
0.0 + 1.0 ε

Añadimos una función de convertir, por defecto

import Base: convert
convert(::Type{Dual}, x::Real) = Dual( (x,zero(x)) );

Por último, tenemos que explicar a Julia que si intenta sumar un real y un dual, hay que “promocionar” el real a dual

import Base: promote_rule
promote_rule(::Type{Dual}, ::Type{<:Number}) = Dual;

Vemos que esto ya nos permite interoperar con Float y Dual

1.0/(1.0+ϵ)
1.0 + -1.0 ε

Cálculo formal de derivadas

Con la construcción que hemos hecho, si f es una función racional \[ g(x + \varepsilon) = g(x) + g'(x) \varepsilon \]

D(g,x) = (g(x + ϵ)).f[2];

Utilizando estos números duales, muchas funciones se pueden derivar de esta forma.

Por ejemplo, tomemos la raíz cuadrada con la formula babilonia

function sqrtBab(x; N=10)
    t = (1 + x)/2
    for i=2:N t=(t + x/t)/2 end
    return t
end

sqrtBab(1.35) - sqrt(1.35)
0.0

Y si aplicamos números duales tenemos que coincide con la derivada simbólica

Dsqrt(x) = .5/√x 
Dsqrt(1.35) - D(sqrtBab,1.35)
-5.551115123125783e-17

Nativos en Julia

Los números duales están implementados en el paquete DualNumbers.jl

Hablaremos más adelante del paquete ForwardDiff.jl

El algoritmo de Euclides

Para números enteros

function gcd(a,b)
    if iszero(b)
        return a
    else
        return gcd(b,mod(a,b))
    end
end

gcd(30,15)
15

Extensión a polinomios

Podemos extender gcd a cualquier lugar con iszero y mod.

Queremos extender la funciones anteriores anterior, de modo que gcd funcione sin cambios

Por esto julia es uno de los lenguajes con más reciclaje de código.

Estructuras

Creamos una estructura para polinomios de coeficientes reales.

Basta con almacenar los coeficientes

struct Poly
    coeff::Vector{Float64} # [a₀, a₁,... , aₙ]

Estructuras

Creamos una estructura para polinomios, que almacena los coeficientes.

struct Poly
    coeff::Vector{Float64} # [a₀, a₁,... , aₙ]
    
    # Constructor strips zeros before leading term
    function Poly(v)
        iNonZero = length(v);
            
        while (iNonZero >= 1) && iszero(v[iNonZero])
            iNonZero = iNonZero - 1;
        end
        new(v[1:iNonZero]);
    end
end

Añadimos un constructor que elimina los ceros

Lo probamos

P = Poly([0,1])
Poly([0.0, 1.0])
P = Poly([0,1,0])
Poly([0.0, 1.0])

Funciones sobre polinomios

function degree(P::Poly)
    degree = length(P.coeff) - 1
end
degree (generic function with 1 method)
function lead(P::Poly)
    if degree(P) < 0 
        return 0
    else 
        return P.coeff[degree(P)+1]
    end
end
lead (generic function with 1 method)
P = Poly([0,1,0])
print("degree(P) = ", degree(P), ", lead(P) = ", lead(P))
degree(P) = 1, lead(P) = 1.0

Haciéndolo bonito

function Base.show(io::IO,P::Poly) 
    if degree(P) < 0
        text = "0"
    else
        text = "$(P.coeff[1])"
        for n = 2:degree(P) + 1
            text = text * " + $(P.coeff[n])" * "x^$(n-1)"
        end
    end
    print(text)
end;
print(P)
0.0 + 1.0x^1

Multiple dispach

Vamos a hacer una función iszero para polinomios

import Base: iszero 

function iszero(P::Poly)
    return iszero(sum(abs.(P.coeff)))
end
iszero (generic function with 16 methods)

Multiple dispatch. Operaciones

import Base: +
function +(a::Poly,b::Poly)
    da = degree(a); db=degree(b)
    return Poly( 
        [a.coeff; zeros(max(db-da,0))] 
        + [b.coeff; zeros(max(da-db,0))]
        )
end
+ (generic function with 266 methods)
P = Poly([0,1,2])
print("P+P = ", P+P)
P+P = 0.0 + 2.0x^1 + 4.0x^2

Multiple dispatch. Operaciones

import Base: *
function *(a::Poly,b::Poly)
    da = degree(a); db = degree(b);
    prodcoeffs = zeros(da+db+1);
    for k=0:(da+db)
        rangei = max(0,k-db):min(k,da)
        prodcoeffs[k+1] = sum(  a.coeff[1   .+ rangei] .*
                                b.coeff[k+1 .- rangei]);
    end
    return Poly(prodcoeffs)
end
* (generic function with 373 methods)
import Base: -
function -(a::Poly,b::Poly)
    return a + Poly([-1])*b
end
- (generic function with 273 methods)
P = Poly([0,1,2])
print("P*P = ", P*P, " P-P = ", P-P)
P*P = 0.0 + 0.0x^1 + 1.0x^2 + 4.0x^3 + 4.0x^4 P-P = 0

Multiple dispatch. mod

import Base: mod
function mod(a::Poly,b::Poly)
    
    if degree(a) < 0 # El polinomio de grado negativo es el 0
        error("¡No dividas por 0!")
    end

    r = a             # En cada paso a = b × q + r
    
    while degree(r) ≥ degree(b) 
        s = Poly([ zeros(degree(r)-degree(b)) ; lead(r)/lead(b) ])
        r = r - s*b
    end
    
    return r
end
mod (generic function with 25 methods)
Q = Poly([-3,1])
P = Q*Q*Q + Poly([2])

print("mod(P,Q) = ",mod(P, Q))
mod(P,Q) = 2.0

Grand finale

Lanzamos la función gcd ¡que nunca oyó hablar de polinomios! Y obtenemos

Q = Poly([-3,1])
P = Q*Q*Q

print("gcd(P,Q)=",gcd(P,Q))
gcd(P,Q)=-3.0 + 1.0x^1