Efficient implementation of covariance updates

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

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

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