4
$\begingroup$

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. [1]

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.

[1] https://en.wikipedia.org/wiki/Centering_matrix

$\endgroup$
5
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ "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! $\endgroup$ – mweylandt Aug 30 '18 at 18:09
  • $\begingroup$ If you'd tell us more about the size and sparsity of your $A$ matrices we might be able to provide more targeted advice. $\endgroup$ – Brian Borchers Aug 30 '18 at 19:40
  • $\begingroup$ 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). $\endgroup$ – mweylandt Aug 31 '18 at 0:16
  • $\begingroup$ I'd look at LSQR as an iterative solver. $\endgroup$ – Brian Borchers Aug 31 '18 at 0:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.