Efficient implementation of covariance updates

How do I implement the following efficiently? [X \in R^{N \times p}, y, \tilde{\beta} \in R^{p}, p \gg N ]

[f(X, y, w, \tilde{\beta}) = \sum_{k:|\tilde{\beta}_k | greater 0 }} (x_j, x_k) \tilde{\beta}_k ]

for i in n_iterations
    for j in p
        f(X,y,w, beta)
    end
end