# Method to Efficiently Solve “Centered” Least Squares without centering “A”

Suppose I want to solve

$$\text{arg min}_x \frac{1}{2}\|\tilde{A}x - b\|_2^2 + \frac{1}{2}\|x - c\|_2^2$$

where $A$ is a wide sparse matrix and $\tilde{A} = A C_n = A (I - \mathbf{1}\, \mathbf{1}^T/n)$ is the "centered" version of $A$, typically non-sparse. 

If I am willing to materialize $\tilde{A}$, I know several efficient methods to solve this problem, but I'd like to avoid doing so if possible. ($\tilde{A}$ is small enough to fit on a single machine, but large enough that a sparsity-aware method may be worthwhile.)

Is there any literature on this problem? (It arises fairly often in statistics, where centering $A$ makes the intercept of a linear regression independent of the other coefficients.)

If it matters, I'll actually need to solve this problem repeatedly for different values of $c$, so I am willing to "pay for" a more expensive factorization and reuse it if possible.

How would you solve the problem if you didn't need to do the centering?

Since $A$ is large and sparse, you'd probably pick an iterative method such as CGNE which depends on being able to perform matrix-vector multiplies $Ax$ and $A^{T}y$. It turns out that you can still use the same iterative method for the centered version of the problem since matrix-vector multiplications with $\tilde{A}$ aren't really slower than matrix-vector multiplications with $A$.

$\tilde{A}x=A(I- 1\; 1^{T}/n)x=Ax-(A1/n)(1^{T}x)$

and

$\tilde{A}^{T}y=(A(I-1\;1^{T}/n)^{T}y=A^{T}y-1(1^{T}A^{T}/n)y$.

In these two expressions, you'd only have to compute $A1/n$ and its transpose $1^{T}A^{T}/n$ once. The rest of the work takes trivial $O(n)$ time.

If you were going to use a QR factorization of $A$, you'd have to deal with the fact that the Q matrix is typically fully dense even though $A$ is sparse. You could use a rank-one update procedure (qrupdate in MATLAB) to update the factorization of $A$ to a factorization of $\tilde{A}$, but this would probably be no faster than centering $A$ and then finding its QR factorization.

• "How would you solve the problem if you didn't need to do the centering?" Not sure - I'm adding support for sparse $A$ to software which currently only supports dense $A$. Centering doesn't hurt anything in the dense case, so I can use standard methods. (In my case, a combination of the Woodbury identity and pre-computing the Cholesky factorization of $I + \tilde{A}\tilde{A}^T$ and using it repeatedly.) Adopting an iterative method in the way you describe makes sense - I'll run some benchmarks and report back in case anyone else stumbles upon this. Thanks! – mweylandt Aug 30 '18 at 18:09
• If you'd tell us more about the size and sparsity of your $A$ matrices we might be able to provide more targeted advice. – Brian Borchers Aug 30 '18 at 19:40
• Size is typically on the order of 1e2-1e4 by 5e2-5e5, with the smaller end of that scale being a bit more typical. Definitely not a scale where a dense approach is impossible, but I'd prefer to do something smart if possible since I'll have desktop users. Sparsity pattern is pretty random: this is statistical software, where $A$ (or $X$ in stat-speak) comes from the data; an example would be a "bag-of-words" model where $A_{ij}$ is either an indicator (whether word $j$ appeared in document $i$) or a not-too-large count. No idea about "typical" spectra or anything like that (but I could look). – mweylandt Aug 31 '18 at 0:16
• I'd look at LSQR as an iterative solver. – Brian Borchers Aug 31 '18 at 0:49