I have an issue with the datafit update step, because the gradient of 0.5 || Xw - y + a 1_n|| with respect to a is n a - (Xw - y)^t 1_n, so to compute it one needs y.sum()

@shz9 do people usually fit an intercept with such models ? How do they do it ?

EDIT: ok I made it work with the user passing y_mean. If not, I don't think it's possible to compute an intercept.

We usually assume $y$ is centered. So, you could go with the implementation with y_mean, but set it to 0 by default?

I just added a way to load the design matrix in a sparse scipy format from magenpy.
Would show us how to load X.T y from magenpy @shz9 ?

Sure, here's a snippet that you can easily integrate for the purposes of testing:

import magenpy as mgp

# Initialize simulator object:
g_sim = mgp.PhenotypeSimulator(mgp.tgp_eur_data_path(), h2=0.5, pi=0.1)
# Simulate outcome (phenotype):
g_sim.simulate()

# Access X:
X_csr = g_sim.genotype[22].to_csr()
# Convert to csc:
X_csc = X_csr.tocsc()
# Access y:
y = g_sim.to_phenotype_table()['phenotype'].values

# Sample size:
n = X_csr.shape[0]

# To get X'y (you may need to standardize both X and y before this step):

Xty = X_csr.T.dot(y)

# For a more general case that will be useful for GWAS settings:

# (1) Perform GWAS:
g_sim.perform_gwas()

# (2) Extract marginal betas (Xty, assumes both x and y are standardized column-wise):
g_sim.sumstats_table[22].standardized_marginal_beta

Would it be possible to have interface take $X^\top X$ instead of $X$ in some instances?

Thanks a lot for the code snippet!

Would it be possible to have the interface take X.T X instead of X in some instances?

This should take more time because we will have to implement a new solver.
We will start with making the code work for X, and then X.T X.

OK, thanks. By the way, I found some interesting statistics literature on losses of the form that we're working on here, going back to Loh and Wainwright (2012):

• High-dimensional regression with noisy and missing data: Provable guarantees with nonconvexity
• CoCoLasso for high-dimensional error-in-variables regression
• An Investigation of Methods for Handling Missing Data with Penalized Regression

These papers were focused on the regime of missing/noisy data in $X$ and they proposed losses similar to the ones we're discussing. So, the work here could have wide applicability if we manage to provide efficient and robust solvers for these losses.

This literature also provides some important caveats and considerations to keep in mind, for example, the fact that $X^\top X$ has to be positive (semi-)definite.