2
$\begingroup$

Is there a permutation used in sparse QR factorizations that better locates small elements on the diagonal of R to the end of the diagonal? As an example, consider the following snippet of MATLAB code:

A = [1  1  1; ...
    -1 -1  1; ...
     0  0 -1; ...
     0  0  1; ...
     -1 0  0; ...
     0 -1  0];
s = [1e-8;1e-8;1e-8;1;1e-8;1];
B = sparse([A diag(s)]);
[Q,R,p] = qr(sparse(B)');
pp = [5,6,1,2,4,3];
[QQ,RR] = qr(sparse(B(pp,:)'));

This matrix is similar to what may occur in an optimization problem with small slack variables. If we look at the absolute value of the diagonal elements of R, we see a small element in position 5:

Diagonal elements of R

I would prefer this element to be at the end of the diagonal of R. Using the permutation pp, this can be achieved and the the diagonal of RR is seen in this image:

enter image description here

In this case, the sparsity between R and RR is the same. Ideally, I'd like a diagonal sorted by absolute value similar to what is achieved in a dense QR factorization. In that case, it's easy to get a rough idea of the effective rank of the original matrix and truncate the matrix R effectively for the applications that I care about. At the same time, I understand that this would dramatically impact the sparsity of R and for large matrices this is not tractable. In truth, I don't really need a fully sorted diagonal, but some kind of permutation that mostly puts the small elements at the end. At the moment, I have some very large matrices and may get some very small elements in the first few dozen diagonal elements, which is undesirable.

As some additional information, if it helps, I understand that MATLAB uses spqr as its underlying QR factorization algorithm and that the parameter spqrtol can be set in spparms to adjust the drop tolerance for small elements. Alternatively, spqr could just be used directly. In this case, I'd contend the issue is how to effectively set this tolerance if the issue is small diagonal elements compared to other diagonal elements.

Thanks for any insights.

$\endgroup$
5
  • $\begingroup$ The diagonal is sorted by absolute value in a dense QR? $\endgroup$ Jul 4 at 19:37
  • $\begingroup$ In MATLAB's case, yes. Per their documentation, Permutation information, returned as a matrix or vector. The shape of P depends on the value of outputForm. Also, qr selects P to satisfy different criteria depending on whether the first input matrix is full or sparse: Full — qr selects P so that abs(diag(R)) is decreasing. Sparse — qr selects P to reduce fill-in in R. Octave does the same, If the matrix A is full, the permuted QR factorization [Q, R, P] = qr (A) forms the QR factorization such that the diagonal entries of R are decreasing in magnitude order. $\endgroup$
    – wyer33
    Jul 4 at 19:56
  • $\begingroup$ As a follow-up, I believe Octave uses geqp3 to compute the dense factorization and likely MATLAB does as well. It's not clear to me from the routine's documentation as to the sorting method, but Octave appears to take the result more or less directly and the diagonal elements are sorted. $\endgroup$
    – wyer33
    Jul 4 at 20:03
  • $\begingroup$ That's a different algorithm though, QRP or rank-revealing QR, which is called only if you request the three-output form explicitly. If you just write [Q, R] = qr(A) it does regular QR, and the entries on diag(R) are not necessariy sorted. $\endgroup$ Jul 4 at 20:25
  • $\begingroup$ That's true and I'm sorry for the confusion. The above question refers to rank-revealing QR factorizations of which spqr provides a sparse and geqp3 provides a dense. I'm interested in whether there's a permutation for sparse matrices that better moves the small diagonal elements to the end of the diagonal of R and still provides reasonable sparsity. This can be done with a rank-revealing QR factorization on dense matrices where density is not an issue. I don't require a fully sorted diagonal, just something to get the small elements farther down the diagonal. $\endgroup$
    – wyer33
    Jul 4 at 20:39

0

Your Answer

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

Browse other questions tagged or ask your own question.