Least Squares with Dense-Block Diagonal Structure

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.

• You've used $p$ for two things (the left hand side of the arg min and a dimension of $V$.) Does $V$ have more rows than columns or vice versa? – Brian Borchers Jun 29 '18 at 4:37
• Sorry for the ambiguities, I fixed those. $V$ has more rows than columns. – Spenser Jun 29 '18 at 4:50
• How big are $m$ and $N$? – Brian Borchers Jun 29 '18 at 4:57
• By the way, what do you mean by "block diagonal" $V$? Does $V$ have non-zero square blocks on the diagonal in its first $m$ rows and $0$'s in all rows $m+1$ through $N$? If so, you should realize that only the first $m$ rows of $J$ are relevant? – Brian Borchers Jun 29 '18 at 23:06
• $N$ is much larger than $m$. And the block-diagonal structure is such that all rows of $V$ have nonzero elements. What I mean by block diagonal is that it can be arranged so that $V$ is a diagonal block matrix, if the blocks in the block matrix are allowed to have arbitrary (not the same) dimensions. – Spenser Jun 30 '18 at 5:56

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.

• Great, thanks! Yes the product $JV$ can be formed while exploiting the diagonal structure, but my primary concern is the actual least-squares solve (i.e. factorization) step itself and whether that can be accelerated by exploiting this structure using something like an iterative solve (which you addressed). It seems that you're saying that to your knowledge the usual QR using the dense $JV$ product may be as good as it gets. – Spenser Jul 1 '18 at 21:05
• Convention iterative methods for least squares problems like LSQR work well for sparse problems but can’t compete on dense problems. There is some relevant research on randomized algorithms for approximate solutions to dense overdetermined least squares problems- see The Blendenpick algorithm for example. However, your problem is small enough that this unlikely to be terribly helpful. – Brian Borchers Jul 2 '18 at 13:57