Computing $\sum_i f_i \mathbf{x}_i \mathbf{x}_i^T$ which is the form of Hessian for Poisson regression

In [1]:

N = int(1e5);
d = 100;
f = randn(N);
x = randn(N, d);
H = zeros(d, d);

In [2]:

function hess1(x, f)
    N, d = size(x);
    temp = zeros(N, d);
    @simd for kk = 1:N
        @inbounds temp[kk, :] = f[kk] * x[kk, :];
    end
    H = x' * temp;
end

Out[2]:

hess1 (generic function with 1 method)

In [3]:

function hess2(x, f)
    N, d = size(x);
    H2 = zeros(d,d);
    @simd for k = 1:N
        @inbounds H2 += f[k] * x[k, :]' * x[k, :];
    end
    return H2
end

Out[3]:

hess2 (generic function with 1 method)

In [4]:

function hess3(x, f)
    N, d = size(x);
    H3 = zeros(d,d);
    for k = 1:N
        for k1 = 1:d
            @simd for k2 = 1:d
                @inbounds H3[k1, k2] += x[k, k1] * x[k, k2] * f[k];
            end
        end
    end
    return H3
end

Out[4]:

hess3 (generic function with 1 method)

In [13]:

@time H1 = hess1(x, f);
@time H2 = hess2(x, f);
@time H3 = hess3(x, f);

elapsed time: 0.776116469 seconds (262480224 bytes allocated, 26.49% gc time)
elapsed time: 30.496472345 seconds (16385442496 bytes allocated, 56.07% gc time)
elapsed time: 2.769934563 seconds (80128 bytes allocated)

In [9]:

norm(H1 - H3)

Out[9]:

6.830690109677714e-11