Least Squares with Dense-Block Diagonal Structure

Question

I need to solve a least squares problem that takes the following form:

$$p = \arg \min_{x}\Vert J V x - y \Vert_2, $$

where $J \in \mathbb{R}^{N \times N}$ is a general dense matrix, and $V \in \mathbb{R}^{N \times m}$ is a block-diagonal matrix with $m << N$ (that is, this is a typical overdetermined least-squares problem).

I'm wondering if there is a way to take advantage of that block-diagonal matrix. When I form the product $JV$, now it is a dense matrix that I can use a dense least-squares solver to solve, but it feels like there should be a better solution. Does anyone have an idea to do this faster?

Edits for additional details: Typical values of $V$ for my current problem are about 100,000 rows and 100 columns, though I'd like to scale that up to many more rows and probably a few more columns. The block diagonal structure means that this matrix could be written as a block matrix that is diagonal, but the blocks are not necessarily square or the same size.

Answer 1

If $N$ is on the order of 100,000 and $m$ is on the order of $100$, Then $J$ requires about 80 gigabytes to store in double precision and $V$ requires a trivial amount of storage. The product $M=JV$ is of size $N$ by $m$ and would be fully dense, requiring about 80 megabytes to store. You should have no trouble storing $V$ or $M$, but you may have to keep $J$ out on disk unless you've got a lot of RAM available.

Once you have $M=JV$ as a matrix, solving the $100,000$ by $100$ linear least squares problem using a QR factorization is quite easy (takes about 1/2 of a second using LAPACK to do the QR factorization on my desktop machine.) Note that you want to use a compact or "Q-less" form of the QR factorization to avoid forming an $N$ by $N$ dense $Q$ matrix.

You can take advantage of the structure of $V$ in multiplying out $M=JV$. If $V$ has $r$ rows and columns of blocks, with only the diagonal blocks being nonzero, then you can partition $J$ in corresponding fashion and do the block matrix multiplication as

$M=JV=\left[ \begin{array}{} J_{1,1} V_{1,1} & J_{1,2} V_{2,2} & \ldots & J_{1,r}V_{r,r} \\ J_{2,1} V_{1,1} & J_{2,2} V_{2,2} & \ldots & J_{2,r}V_{r,r} \\ \vdots & \vdots & \vdots & \vdots \\ J_{r,1} V_{1,1} & J_{r,2} V_{2,2} & \ldots & J_{r,r} V_{r,r} \end{array} \right] $

This can be done one block column at a time, and if $J$ is stored in a file would require only one pass through the file.

Iterative methods are not a good idea for this problem. Although you could do multiplications of $V$ and $V^{T}$ times a vector easily, multiplications involving $J$ would be very expensive. If you chose to just multiply $M=JV$ once, then you're faced with a dense 100,000 by 100 linear least squares problem which is faster to solve by QR factorization than by iterative methods.