sparse QRCP example: unusual LU-based co-range stabilization

[rileyjmurray]

RandBLAS' example for truncated QRCP of sparse matrices employs power iteration $W \mapsto W A$ on wide matrices $W$. That power iteration can be stabilized by Householder LQ, sketch-orthogonal Cholesky LQ, or an LU-based method. The implementation of the LU-based method (which I'm responsible for, below ...) is odd.

RandBLAS/examples/sparse-low-rank-approx/qrcp_matrixmarket.cc

Lines 177 to 200 in 71b1b9d

	int lu_row_stabilize(int64_t m, int64_t n, T* mat, int64_t* piv_work) {
	randblas_require(m < n);
	for (int64_t i = 0; i < m; ++i)
	piv_work[i] = 0;
	lapack::getrf(m, n, mat, m, piv_work);
	// above: the permutation applied to the rows of mat doesn't matter in our context.
	// below: Need to zero-out the strict lower triangle of mat and scale each row.
	T tol = std::numeric_limits<T>::epsilon()*10;
	bool nonzero_diag_U = true;
	for (int64_t j = 0; (j < m-1) & nonzero_diag_U; ++j) {
	nonzero_diag_U = abs(mat[j + j*m]) > tol;
	for (int64_t i = j + 1; i < m; ++i) {
	mat[i + j*m] = 0.0;
	}
	}
	if (!nonzero_diag_U) {
	throw std::runtime_error("LU stabilization failed. Matrix has been overwritten, so we cannot recover.");
	}
	for (int64_t i = 0; i < m; ++i) {
	T scale = 1.0 / mat[i + i*m];
	blas::scal(n, scale, mat + i, m);
	}
	return 0;
	}

The method starts by decomposing $P W = LU$ for a permutation $P$ and unit lower-triangular $U$. This is immediately a weird thing to do because the pivot decisions from running LU on an $m \times n$ matrix only depend on the first $\min\{m,n\}$ columns of that matrix. Then it implicitly forms a diagonal matrix $D$ where $D_{ii} = U_{ii}^{-1}$, raising an error if $U_{ii} =0$, and replaces $W \leftarrow D^{-1} U$.

A proper way of doing LU-based co-range stabilization is to decompose $P W^T = L U$, so that $W = U^T L^T P$, then returning $W \leftarrow L^T P$; see here for an implementation. This procedure should never raise an error. The only thing of note is that if $W$'s rows are linearly dependent on entry then they'll be made linearly independent on exit, which means the co-range stabilization can spuriously increase the dimension of the subspace.