Оптимизировать сильно векторизованный код в Юлии?

Вот код, который все еще векторизован, и он не использует многопоточность. На моем компьютере он более чем в два раза быстрее оригинала, но кое-что еще можно сделать здесь.

Кстати, я серьезно сомневаюсь, что код Matlab особенно оптимизирован, так как происходят некоторые очень расточительные операции - ненужное распределение, ненужные операции (sum(w.*1) действительно плохо, умножение массива на 1 и выделение дополнительного массива в процессе это ужасно расточительно :) Кроме того, вам не нужно выделять какие-либо векторы xx_0_0, xx_1_0 в Matlab, вы можете просто использовать трансляцию, как в Джулии.

Во всяком случае, вот моя первая попытка:

function fit_xlin2(x, y, w)
    regularization = 1.0e-5

    sumwx = (w' * x) + regularization
    sumwy = (w' * y)
    sumwxx = sum(a[1]*a[2]^2 for a in zip(w, x))
    sumwxy = sum(prod, zip(w, x, y))

    wx = w .* x
    xxk_0_0 = sum(w) .- w
    xxk_1_0 = sumwx .- wx
    xxk_1_1 = sumwxx .- wx .* x
    xyk_0 = sumwy .- w .* y
    xyk_1 = sumwxy .- wx .* y

    det = xxk_0_0 .* xxk_1_1 .- xxk_1_0 .* xxk_1_0
    c0  = (xxk_1_1 .* xyk_0  .- xxk_1_0 .* xyk_1)./det
    c1  = (-xxk_1_0 .* xyk_0 .+ xxk_0_0 .* xyk_1)./det

    return c0 .+ c1 .* x
end

Редактировать: Вы можете получить достаточное ускорение, если де-векторизовать основной цикл. Этот код примерно в 17 раз быстрее исходного кода Julia, все еще однопоточный и достаточно читаемый:

function fit_xlin_loop(x, y, w)
    if !(size(x) == size(y) == size(w))
        error("Input vectors must have the same size.")
    end

    regularization = 1.0e-5

    sumw = sum(w)
    sumwx = (w' * x) + regularization
    sumwy = (w' * y)
    sumwxx = sum(a[1]*a[2]^2 for a in zip(w, x))
    sumwxy = sum(prod, zip(w, x, y))

    y_est = similar(x)
    @inbounds for i in eachindex(y_est)
        wx = w[i] * x[i]
        xxk_0_0 = sumw - w[i]
        xxk_1_0 = sumwx - wx
        xxk_1_1 = sumwxx - wx * x[i]
        xyk_0 = sumwy - w[i] * y[i]
        xyk_1 = sumwxy - wx * y[i]

        det = xxk_0_0 * xxk_1_1 - xxk_1_0 * xxk_1_0
        c0  = (xxk_1_1 * xyk_0 - xxk_1_0 * xyk_1) / det
        c1  = (-xxk_1_0 * xyk_0 + xxk_0_0 * xyk_1) / det

        y_est[i] = c0 + c1 * x[i]
    end
    return y_est
end