Meu fluxo de trabalho em Julia: Legibilidade VS. Desempenho

Olá!

Após o meu último post, recebi muitas perguntas como: “O que é Julia? Por que você está usando Julia? Vale a pena?”

Bem, vou tentar responder as duas últimas perguntas mostrando o meu fluxo de trabalho: como eu codifico e que tipo de ganhos eu obtenho (para a primeira pergunta eu vou indicar o website oficial 🙂 ).

Requisitos necessários

É necessário um pouco de conhecimento de C e Julia para que você possa entender totalmente este post.

Meu caso de uso

Muitas das minhas atividades no INPE exigem análises que só podem ser realizadas por algoritmos de otimização ou simulações de Monte Carlo. Por exemplo, durante os estudos da Pré-Fase A no Centro de Projeto Integrado de Missões Espaciais (CPRIME), às vezes utilizamos algoritmos evolutivos para obter as melhores soluções para o projeto de um satélite. Também é comum utilizar simulações de Monte Carlo para encontrar os piores cenários possíveis, por exemplo, o estado inicial do satélite (atitude e velocidade angular) que produz o maior tempo para adquirir o apontamento inicial para o Sol. Assim, alguns algoritmos que codificamos devem ser executados milhares, até milhões, de vezes. Então, o desempenho é algo que realmente importa para mim.

Por outro lado, somos principalmente engenheiros que estão acostumados com o MATLAB. Os algoritmos são escritos para resolver equações e simular sistemas dinâmicos. Portanto, quanto mais próximo o algoritmo estiver das equações matemáticas, melhor.

A seguir, descreverei meu fluxo de trabalho: como eu escrevo o código em Julia desde a validação inicial até a versão que deve ser utilizada dentro de um otimizador. Eu tentei muitas abordagens ao longo dos anos, mas esta foi a que forneceu o menor tempo para obter os resultados.

Exemplo numérico

Para ilustrar, vamos supor que queremos um algoritmo para executar a atualização de covariância do filtro de Kalman de um sistema dinâmico. Não dê muita importância a essas equações, apenas selecionei um cenário aleatório para ilustrar melhor o que estamos fazendo.

Inicialização

\(\begin{aligned}\mathbf{F}_{k-1} &= \exp\left[\begin{array}{c c c} \mathbf{w}^{\times} & -\mathbf{I}_3 & \mathbf{0}_{3 \times 12} \\ \mathbf{0}_{3 \times 3} & \lambda \cdot \mathbf{I}_3 & \mathbf{0}_{3 \times 12} \\ \mathbf{0}_{12 \times 3} & \mathbf{0}_{12 \times 3} & \mathbf{0}_{12 \times 12 } \end{array} \right] \\ \mathbf{H}_k &= \left[ \begin{array}{c c c c c c} \mathbf{s}^{\times} & \mathbf{0}_{3 \times 3} & \mathbf{s}^{\times} & \mathbf{I}_3 & \mathbf{0}_{3 \times 3} & \mathbf{0}_{3 \times 3} \\ \mathbf{B}^{\times} & \mathbf{0}_{3 \times 3} & \mathbf{0}_{3 \times 3} & \mathbf{0}_{3 \times 3} & \mathbf{B}^{\times} & \mathbf{I}_3 \end{array} \right] \\ \mathbf{P}_p &= \mathbf{I}_{18} \\ \mathbf{P}_u &= \mathbf{I}_{18} \end{aligned}\)

Loop de simulação

Para \(k = 2 \) até \(60000\), faça:

\(\begin{aligned} \mathbf{P}_p &= \mathbf{F}_{k-1} \mathbf{P}_p \mathbf{F}_{k-1}^{T} + \mathbf{Q} \\ \mathbf{K}_k &= \mathbf{P}_p\mathbf{H}_k^T(\mathbf{R} + \mathbf{H}_{k} \mathbf{P}_p \mathbf{H}_{k}^{T}) \\ \mathbf{P}_u &= (\mathbf{I}_{18} – \mathbf{K}_k\mathbf{H}_k)\mathbf{P}_p(\mathbf{I}_{18}-\mathbf{K}_k\mathbf{H}_k)^T + \mathbf{K}_k\mathbf{R}\mathbf{K}_k^T \end{aligned} \)

Toda a vez que \(\mathbf{P}_u \) for modificada, guarde \(\mbox{trace}(\mathbf{P}_u) \).

Nota

\(\mathbf{w}^{\times} \), \(\mathbf{s}^{\times} \) e \(\mathbf{B}^{\times} \) são matrizes constantes, \(\mathbf{0}_{n \times m} \) é uma matriz \(n \times m\) preenchida com 0s, e \(\mathbf{I}_n\) é uma matrix identidade \(n \times n\).

Primeiro: escreva o código rápido, o mais legível possível

A primeira coisa é codificar uma versão inicial do algoritmo que você quer, sem se preocupar muito com o desempenho. Nesta fase, a legibilidade é fundamental, porque facilitará a depuração.

Nossa primeira versão do código será:

using LinearAlgebra

# Constants
# ==============================================================================

const O3x3   = zeros(3,3)
const O3x12  = zeros(3,12)
const O12x18 = zeros(12,18)

function kalman_filter()

    # Constants
    # ==========================================================================

    λ = 1/100
    Q = 1e-20I
    R = 1e-2I

    # Dynamic model
    # ==========================================================================

    wx = Float64[ 0 -1  0;
                  1  0  0;
                  0  0  0;]

    Ak_1 = [  wx    -I   O3x12;
             O3x3  λ*I   O3x12;
                  O12x18      ;]

    Fk_1 = exp(Ak_1)

    # Measurement model
    # ==========================================================================

    sx = Float64[ 0  0  0;
                  0  0 -1;
                  0  1  0;]

    Bx = Float64[ 0  0 -1;
                  0  0  0;
                  1  0  0;]

    Hk = [sx O3x3   sx    I O3x3 O3x3;
          Bx O3x3 O3x3 O3x3   Bx    I;]

    # Kalman filter initialization
    # ==========================================================================

    Pu = Matrix{Float64}(I,18,18)

    # Simulation
    # ==========================================================================

    result = Vector{Float64}(undef, 60000)

    result[1] = tr(Pu)

    for k = 2:60000
        Pp = Fk_1*Pu*Fk_1' + Q
        K  = Pp*Hk'*pinv(R + Hk*Pp*Hk')
        Pu = (I - K*Hk)*Pp*(I - K*Hk)' + K*R*K'

        result[k] = tr(Pu)
    end

    return result
end

É possível ver quão perto do equacionamento matemático esse código está? Mesmo os símbolos como λ podem ser utilizados em Julia, que suporta qualquer caractere UTF-8 no nome das variáveis. Agora, imagine algo com centenas de equações. É mais fácil validar o código e encontrar problemas quando você escrever usando uma linguagem computacional parecida com a linguagem matemática.

Informação

O I em Julia não é uma matriz identidade 3 x 3 como nas equações acima. Ele é chamado de operador UniformScalling. Em nosso código, ele se transforma em uma matriz identidade com a dimensão necessária para a operação fazer sentido, isto é 3 x 3.

Nesta primeira versão do código, temos duas tarefas: 1) ver se funciona corretamente; e 2) ver se o desempenho é bom o suficiente para as nossas necessidades. Sabemos que após 60.000 iterações, o traço da matriz de covariância \(P_u\) é aproximadamente 7,00089. Em Julia, podemos executar o código no REPL e ver se isso é verdade:

julia> result = kalman_filter();

julia> result[end]
7.000894875771719

Legal 🎊🎉! Pelo menos o código faz o que é necessário. O próximo passo é medir o desempenho e verificar se é o suficiente. Isso é muito importante: Julia é uma linguagem muito rápida para mim não porque é mais rápida que C, mas porque o tempo entre a escrita do código inicial e o resultado final é muito menor. Às vezes você simplesmente não precisa de um código com grande desempenho e quer codificar o mais rápido possível, Julia pode fazer isso conforme foi mostrado. Por outro lado, algumas vezes você realmente precisa de velocidade, Julia pode fazer isso também como será mostrado a seguir.

A maneira correta de avaliar o desempenho de um código em Julia é colocá-lo dentro de uma função (já fizemos isso) e usar o excelente pacote BenchmarkTools.jl:

julia> using BenchmarkTools

julia> @btime kalman_filter()
  1.296 s (2640271 allocations: 2.83 GiB)

Agora sabemos que o tempo médio de execução da função (após a etapa de compilação) é de 1,296s e que foi alocado 2,83 GiB de memória. Neste momento, eu me pergunto: “Essa velocidade é suficiente para o que eu preciso fazer?” Por exemplo, se eu precisar executar esse código apenas uma vez, posso conviver com seu estado atual. No entanto, se eu precisar que esse código seja utilizado dentro de um otimizador sendo executado 1.000.000 de vezes, eu realmente desejaria otimizá-lo.

Segundo: Modifique seu código para melhorar o desempenho, se necessário

Há algumas coisas que você precisa tomar cuidado para aumentar o desempenho de um código em Julia. A documentação é muito boa e eu realmente sugiro que você a leia. Para começar, você deve se preocupar com três coisas:

  • Estabilidade do tipo: Julia deve saber o tipo das variáveis no código a partir apenas do tipo das entradas. Isso pode ser verificado com a macro @code_warntype.
  • Alocações: o processo para alocar memória é lento, você deve realizar o mínimo de alocações possíveis para melhorar o desempenho.
  • Evite cálculos duplicados: os compiladores estão ficando cada vez mais inteligentes e às vezes eles conseguem identificar cálculos duplicados, mas não faz mal ajudá-los!

Podemos ver facilmente que nosso código é tipo-estável (type stable) ao chamar @code_warntype kalman_filter(). Você verá que nenhuma variável pode ter tipos distintos. Eu não vou me aprofundar aqui porque este é um assunto para um post inteiro. Para maiores informações, consulte a documentação.

Existem alguns métodos para verificar quais partes do código estão alocando. O método “oficial” é iniciar Julia com a opção --track-allocation=user, executar a função e, em seguida, analisar o arquivo .mem que será criado. Isso é descrito aqui. No entanto, existe um método que eu, pessoalmente, acho melhor. Ele é baseado no pacote TimerOutputs.jl. Esse pacote fornece uma maneira de analisar apenas linhas específicas. No nosso caso, devemos nos preocupar com o código executado dentro do loop for. O arquivo anterior deve ser modificado da seguinte maneira:

using LinearAlgebra
using TimerOutputs

...

    reset_timer!()
    for k = 2:60000
        @timeit "Pp" Pp = Fk_1*Pu*Fk_1' + Q
        @timeit "K"  K  = Pp*Hk'*pinv(R + Hk*Pp*Hk')
        @timeit "Pu" Pu = (I - K*Hk)*Pp*(I - K*Hk)' + K*R*K'

        result[k] = tr(Pu)
    end

    print_timer()

    return result
end

Agora, quando você executar kalman_filter();, você verá na tela:

 ──────────────────────────────────────────────────────────────────
                           Time                   Allocations
                   ──────────────────────   ───────────────────────
 Tot / % measured:      2.01s / 69.0%           2.86GiB / 99.2%

 Section   ncalls     time   %tot     avg     alloc   %tot      avg
 ──────────────────────────────────────────────────────────────────
 K          60.0k    673ms  48.5%  11.2μs   1.08GiB  38.3%  19.0Ki
 Pu         60.0k    511ms  36.9%  8.52μs   1.29GiB  45.4%  22.5KiB
 Pp         60.0k    203ms  14.7%  3.39μs    472MiB  16.3%  8.06KiB
 ─────────────────────────────────────────────────────────────────

OK, agora temos algo para analisar. Vamos realizar quatro modificações:

  • Essas três variáveis Pp, K, e Pu estão sendo alocadas a cada iteração. Uma coisa que devemos fazer é alocá-las antes do loop e dizer para Julia utilizar o mesmo espaço para armazenar os novos valores. A última ação pode ser feita utilizando o operador .= que executa um broadcasted assignment. Em outras palavras, cada elemento da variável do lado esquerdo será igual ao elemento correspondente da expressão do lado direito, o que não realiza alocações. Também é bom substituir o operador + pelo operador .+ pela mesma razão, uma vez que a soma de matrizes é uma operação executada elemento por elemento.
  • A variável UniformScaling será trocada por matrizes pré-alocadas. Isso também ajudou a reduzir as alocações.
  • Note que (I - K*Hk) está sendo calculado duas vezes dentro do loop. Nós usaremos uma variável auxiliar para evitar isso. É importante que essa variável auxiliar também seja pré-alocada antes do loop pela mesma razão que discutimos no primeiro ponto.
  • Como podemos garantir que todos os elementos das matrizes dentro do for são acessados por índices válidos, então podemos utilizar a macro @inbounds como descrito na documentação. Isso evita operações de verificação de dimensões, melhorando o desempenho.

Dado isto, nosso código agora deve ser:

using LinearAlgebra

# Constants
# ==============================================================================

const O3x3   = zeros(3,3)
const O3x12  = zeros(3,12)
const O12x18 = zeros(12,18)

function kalman_filter()

    # Constants
    # ==========================================================================

    I6  = Matrix{Float64}(I,6,6)
    I18 = Matrix{Float64}(I,18,18)
    λ   = 1/100
    Q   = 1e-20I18
    R   = 1e-2I6

    # Dynamic model
    # ==========================================================================

    wx = Float64[ 0 -1  0;
                  1  0  0;
                  0  0  0;]

    Ak_1 = [  wx     -I   O3x12;
             O3x3   λ*I   O3x12;
                   O12x18     ;]

    Fk_1 = exp(Ak_1)

    # Measurement model
    # ==========================================================================

    sx = Float64[ 0  0  0;
                  0  0 -1;
                  0  1  0;]

    Bx = Float64[ 0  0 -1;
                  0  0  0;
                  1  0  0;]

    Hk = [sx O3x3   sx    I O3x3 O3x3;
          Bx O3x3 O3x3 O3x3   Bx    I;]

    # Kalman filter initialization
    # ==========================================================================

    Pu = Matrix{Float64}(I,18,18)
    Pp = similar(Pu)

    # Simulation
    # ==========================================================================

    result = Vector{Float64}(undef, 60000)

    result[1] = tr(Pu)

    # Auxiliary variables.
    K    = zeros(18,6)
    aux1 = zeros(18,18)

    @inbounds for k = 2:60000
        Pp   .= Fk_1*Pu*Fk_1' .+ Q
        K    .= Pp*Hk'*pinv(R .+ Hk*Pp*Hk')
        aux1 .= I18 .- K*Hk
        Pu   .= aux1*Pp*aux1' .+ K*R*K'

        result[k] = tr(Pu)
    end

    return result
end

O benchmark dessa nova versão é:

julia> @btime kalman_filter();
  1.109 s (2340283 allocations: 2.06 GiB)

Acabamos de reduzir o tempo computacional em 14% com apenas pequenas modificações! Se essa função for utilizada em um otimizador que a executa 1.000.000 de vezes, economizaremos mais de 2 dias.

Terceiro: Indo ao extremo!

Se é realmente necessário reduzir o tempo computacional, então ainda podemos melhorar algumas coisas aqui. No entanto, alerta de spoiler, isso prejudicará a legibilidade. Por isso, aconselho seguir este caminho somente depois que seu código estiver testado e validado e se você realmente precisar de mais velocidade.

Olhando para o último benchmark, ainda vemos que o código está alocando 2GiB. Isso vem das matrizes temporárias que Julia deve alocar para armazenar multiplicações intermediárias. Por exemplo, na multiplicação Fk_1*Pu*Fk_1' .+ Q, é necessário alocar pelo menos um array temporário para armazenar Fk_1*Pu ou Pu*Fk_1'. O compilador não pode (ainda?) adivinhar que isso pode ser alocado antes do loop ao invés de alocá-lo a cada iteração. Portanto, se realmente precisamos reduzir o tempo de computação, devemos reduzir drasticamente essas alocações.

A abordagem é bastante semelhante para a segunda fase: criar matrizes temporárias fora do loop e usá-las para armazenar as operações intermediárias. No entanto, aqui não podemos usar o operador de broadcast de multiplicação .*, pois ele executará uma multiplicação de elementos ao invés de um produto de matrizes. Dessa forma, a solução é usar a seguinte função:

mul!(Y, A, B)

que calcula o produto matricial A*B e armazena na matriz Y. Essa função também retorna a matriz Y para facilitar a concatenação de operações. Finalmente, o código modificado é:

using LinearAlgebra

# Constants
# ==============================================================================

const O3x3   = zeros(3,3)
const O3x12  = zeros(3,12)
const O12x18 = zeros(12,18)

function kalman_filter()

    # Constants
    # ==========================================================================

    I6  = Matrix{Float64}(I,6,6)
    I18 = Matrix{Float64}(I,18,18)
    λ   = 1/100
    Q   = 1e-20I18
    R   = 1e-2I6

    # Dynamic model
    # ==========================================================================

    wx = Float64[ 0 -1  0;
                  1  0  0;
                  0  0  0;]

    Ak_1 = [  wx     -I   O3x12;
             O3x3   λ*I   O3x12;
                   O12x18     ;]

    Fk_1 = exp(Ak_1)

    # Measurement model
    # ==========================================================================

    sx = Float64[ 0  0  0;
                  0  0 -1;
                  0  1  0;]

    Bx = Float64[ 0  0 -1;
                  0  0  0;
                  1  0  0;]

    Hk = [sx O3x3   sx    I O3x3 O3x3;
          Bx O3x3 O3x3 O3x3   Bx    I;]

    # Kalman filter initialization
    # ==========================================================================

    Pu = Matrix{Float64}(I,18,18)
    Pp = similar(Pu)

    # Simulation
    # ==========================================================================

    result = Vector{Float64}(undef, 60000)

    result[1] = tr(Pu)

    # Auxiliary variables.
    K    = zeros(18,6)
    aux1 = zeros(18,18)
    aux2 = zeros(18,18)
    aux3 = zeros(18,6)
    aux4 = zeros(6,6)
    aux5 = zeros(18,18)
    aux6 = zeros(18,18)

    for k = 2:60000
        # Pp = Fk_1*Pu*Fk_1' + Q
        mul!(aux2, mul!(aux1, Fk_1, Pu), Fk_1')
        @. Pp = aux2 + Q

        # K = Pp*Hk'*pinv(R + Hk*Pp*Hk')
        mul!(aux4, Hk, mul!(aux3, Pp, Hk'))
        mul!(K, aux3, pinv(R + aux4))

        # Pu = (I - K*Hk)*Pp*(I - K*Hk)' + K*R*K'
        mul!(aux1, K, Hk)
        @. aux2 = I18 - aux1
        mul!(aux6, mul!(aux5, aux2, Pp), aux2')
        mul!(aux5, mul!(aux3, K, R), K')
        @. Pu = aux6 + aux5

        result[k] = tr(Pu)
    end

    return result
end

Bem … isso é … feio 😅 No entanto, vamos fazer um benchmark:

julia> @btime kalman_filter()
  848.334 ms (1680299 allocations: 915.17 MiB)

Wow! Esse código é 23% mais rápido que a versão anterior e 34% mais rápido que a primeira versão. Comparado com esta última, a nova versão nos pouparia mais de 5 dias se executássemos o algoritmo 1.000.000 de vezes. O trade-off aqui é ter um código menos legível, matematicamente falando.

Informação

Julia tem excelentes capacidades de meta-programação. Pode-se facilmente construir uma macro que irá expandir as multiplicações para usar as variáveis auxiliares, mantendo a legibilidade. No entanto, este é um tópico um pouco avançado para cobrir por agora.

OK Ronan, e quanto ao C ++?

Essa é a pergunta mais frequente que me fazem quando falo de Julia. Naturalmente, o C ++ pode ser mais rápido no tempo de execução. No entanto, como Julia é uma linguagem dinâmica (parece interpretada), é mais fácil codificar e testar os resultados. Essa é uma das razões pelas quais linguagens como MATLAB e Python foram criadas. Assim, se você considerar o tempo total (desenvolvimento, teste e execução), na minha experiência, você pode vencer o C++ usando Julia (a menos que você seja um mágico em C++, o que a maioria dos engenheiros que eu conheço, inclusive eu, não são).

De qualquer forma, vamos apenas escrever o mesmo algoritmo em C++ para ver o quão longe estamos em termos de desempenho. Observe que estas são minhas habilidades em C++. Tenho quase certeza de que esse código pode ser otimizado (alguém pode fazer isso pra mim, por favor?). No entanto, eu não sou um programador experiente em Julia também. Assim, tenho a sensação de que o código em Julia também pode ser otimizado (alguém, também, pode fazer isso pra mim, por favor?).

A versão em C++ utiliza a biblioteca Eigen para as operações matriciais. Caso contrário, será um código muito, muito longo. O código pode ser visto a seguir:

#include <iostream>
#include "Eigen/Core"
#include "Eigen/QR"
#include "unsupported/Eigen/MatrixFunctions"

using namespace Eigen;

int main(int argc, char *argv[])
{
    double l = 1.0/100.0;
    Matrix<double, 18, 18> Q;
    Matrix<double,  6,  6> R;
    Matrix<double, 18, 18> Ak_1;
    Matrix<double, 18, 18> Fk_1;
    Matrix<double,  6, 18> Hk;
    Matrix<double, 18, 18> Pp;
    Matrix<double, 18, 18> Pu;
    Matrix<double, 18,  6> K;

    VectorXd result(60000);

    // Auxiliary variables.
    MatrixXd I6(6,6);
    MatrixXd I18(18,18);

    MatrixXd aux1(6,6);
    MatrixXd aux2(18,18);
    MatrixXd pinv(6,6);

    // Constants
    // =========================================================================

    I6  = MatrixXd::Identity(6,6);
    I18 = MatrixXd::Identity(18,18);
    Q   = I18*1e-20;
    R   = I6*1e-2;

    // Dynamic model
    // =========================================================================

    Ak_1 << 0, -1,  0, -1,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            1,  0,  0,  0, -1,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  l,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  l,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  l, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0,  0,  0,  0,  0,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;

    // WARNING: Matrix exponential is unsupported in Eigen 3.3.7.
    Fk_1 = Ak_1.exp();

    // Measurement model
    // =========================================================================

    Hk << 0,  0,  0,   0, 0, 0,  0,  0,  0,   1, 0, 0,   0,  0,  0,  0, 0, 0,
          0,  0, -1,   0, 0, 0,  0,  0, -1,   0, 1, 0,   0,  0,  0,  0, 0, 0,
          0,  1,  0,   0, 0, 0,  0,  1,  0,   0, 0, 1,   0,  0,  0,  0, 0, 0,
          0,  0, -1,   0, 0, 0,  0,  0,  0,   0, 0, 0,   0,  0, -1,  1, 0, 0,
          0,  0,  0,   0, 0, 0,  0,  0,  0,   0, 0, 0,   0,  0,  0,  0, 1, 0,
          1,  0,  0,   0, 0, 0,  0,  0,  0,   0, 0, 0,   1,  0,  0,  0, 0, 1;

    // Kalman filter initialization
    // =========================================================================

    Pp.setIdentity();
    Pu.setIdentity();

    // Simulation
    // =========================================================================

    result[0] = Pu.trace();

    for(int k = 2; k <= 60000; k++) {
        Pp   = Fk_1*Pu*Fk_1.transpose() + Q;
        aux1 = R + Hk*Pp*Hk.transpose();
        pinv = aux1.completeOrthogonalDecomposition().pseudoInverse();
        K    = Pp*Hk.transpose()*pinv;
        aux2 = I18 - K*Hk;
        Pu   = aux2*Pp*aux2.transpose() + K*R*K.transpose();

        result[k-1] = Pu.trace();
    }

    return 0;
}

Na minha máquina, o tempo médio de execução utilizando -O3 é:

$ time ./a.out
./a.out  0,83s user 0,00s system 99% cpu 0,836 total

que é muito, muito próximo da versão mais rápida do código em Julia.

Informação

Alguns podem argumentar que seria mais rápido se eu usasse -march=native. De fato, é! No entanto, para uma comparação justa, eu teria também que recompilar a imagem de sistema de Julia com a mesma flag, mas, infelizmente, não tenho tempo para fazê-lo agora. Este é um procedimento fácil se você usar o pacote PackageCompiler.jl.

Conclusão

Espero que você tenha entendido como Julia está me ajudando a atingir meus objetivos mais rapidamente. Agora, temos uma linguagem fácil para codificar e testar algoritmos e que, com pequenas adaptações, pode atingir o desempenho de C. Eu não tenho mais necessidade de traduzir todos os algoritmos que precisam de alto desempenho para uma linguagem compilada, o que acaba sendo muito demorado para um engenheiro que não possui conhecimento avançadíssimo de C / C++.