Memory efficiency (3)

[gdalle]

Julia_AR4_Bayesian_Regression/src/NewB.jl

Lines 38 to 50 in 9e5aa17

	function gibbs(Y,X,BETA0,Sigma0,sig2_d,d0,nu0,n_gibbs,burn)
	beta_dist = zeros(5,n_gibbs-burn)
	sigma_dist = zeros(1,n_gibbs-burn)
	for i = 1:n_gibbs
	beta_d = genBeta(X,Y,BETA0,Sigma0,sig2_d)
	sig2_d = genSigma(Y,X,beta_d,nu0,d0)
	if i > burn
	beta_dist[:,i-burn] = beta_d
	sigma_dist[1,i-burn] = sig2_d
	end
	end
	return beta_dist, sigma_dist
	end

Here we also have an example where you could reuse memory easily for the samples, doing something like this (if you rewrite the associated functions in a mutating way, with dot syntax):

function gibbs(Y,X,BETA0,Sigma0,sig2_d,d0,nu0,n_gibbs,burn)
    beta_dist = zeros(5,n_gibbs-burn)
    sigma_dist = zeros(1,n_gibbs-burn)
    beta_d = ... # init with the right size
    sig2_d = ... # init with the right size
    for i = 1:n_gibbs
        genBeta!(beta_d, X, Y, BETA0, Sigma0, sig2_d)
        genSigma!(sigma_d, Y, X, beta_d, nu0, d0)
        if i > burn
            beta_dist[:,i-burn] .= beta_d
            sigma_dist[1,i-burn] .= sig2_d
        end
    end
    return beta_dist, sigma_dist
end

[borisblagov]

I don't know if I am doing something wrong but I expected the largest gains optimizing here by reducing allocations of the vectors (creating beta_d and sig2_d every time). That doesn't seem to be the case.

First I created the mutating genBeta! by simply adding the dot and beta_d as input

function genBeta!(beta_d,X,Y,Beta_prior,Sigma_prior,sig2_d)
    invSig = Sigma_prior^-1
    V = (invSig + sig2_d^(-1)*(X'*X))^-1
    C = cholesky(Hermitian(V)) 
    Beta1 =  V*(invSig*Beta_prior + sig2_d^(-1)*X'*Y)
    beta_d .= Beta1 .+ C.L*randn(5,1)
end

The new genBeta! has one allocation less than non-mutating function, which is what I expected.

julia> @btime genBeta($X,$Y,$BETA0,$Sigma0,$sig2_d,$beta_d);
  3.688 μs (18 allocations: 4.98 KiB)

julia> @btime genBeta!($beta_d,$X,$Y,$BETA0,$Sigma0,$sig2_d);
  3.663 μs (17 allocations: 4.89 KiB)

Changing the gibbs syntax to

function gibbs(Y,X,BETA0,Sigma0,sig2_d,d0,nu0,n_gibbs,burn)
    beta_d = similar(BETA0)
    beta_dist = zeros(5,n_gibbs-burn)
    sigma_dist = zeros(1,n_gibbs-burn)
    for i = 1:n_gibbs
        genBeta!(beta_d,X,Y,BETA0,Sigma0,sig2_d)
        sig2_d = genSigma(Y,X,beta_d,nu0,d0)
        if i > burn
            beta_dist[:,i-burn] .= beta_d
            sigma_dist[1,i-burn] = sig2_d
        end
    end
    return beta_dist, sigma_dist
end

does not markedly improve the allocations (what I expected), nor the time:

julia> @btime gibbs(Y,X,BETA0,Sigma0,sig2_d,d0,nu0,n_gibbs,burn);
  32.703 ms (140006 allocations: 50.75 MiB)

julia> @btime gibbs_old(Y,X,BETA0,Sigma0,sig2_d,d0,nu0,n_gibbs,burn);
  32.858 ms (147006 allocations: 51.39 MiB)

Also also struggled to define the mutating genSigma! function until I realized that scalars cannot be mutated and therefore indexed into. None of the following worked

sig2_d .= 1/only(sig2_inv)
sig2_d[1,1] = 1/only(sig2_inv)
sig2_d[] = 1/only(sig2_inv)
sig2_d[1] = 1/only(sig2_inv)

so I would have to define a 1-element array and change sig2_d accordingly.

[gdalle]

I don't know if I am doing something wrong but I expected the largest gains optimizing here by reducing allocations of the vectors (creating beta_d and sig2_d every time). That doesn't seem to be the case.

See #4

Also also struggled to define the mutating genSigma! function until I realized that scalars cannot be mutated and therefore indexed into. None of the following worked

My bad, I suggested defining it without remembering that it was a scalar ^^ Glad you found out the hard way

[borisblagov]

Yes, diving into the profiler my next step. I rewrote genSigma to mutate a vector with one element and nothing changed much.

There are also other computations that can be avoided, e.g. X'*X and X'*Y. These of course are not langiage specific and I have them in as this was a code meant for teaching (the matlab version) and hence prioritizes clarity vs speed.