Pivoted Cholesky: Why is rank a private field/property?

Specific Domains Numerics
1

Hello! I am posting for the first time so please forgive me if this is naive.

I am using the cholesky function with pivoting to (1) factorize a positive semidefinite matrix, (2) get its rank, and (3) do some follow-up computations. When I ran propertynames on the CholeskyPivoted object, it did not show rank as a property. Eventually I figured out that rank can still be accessed by doing Achol.rank. See examples:

n = 20; r = 5 # dimensions of matrix

X = reshape(randn(n * r), n,r)
A = X*X' # rank deficient symmetric matrix

Achol = cholesky(A, RowMaximum(); check = false, tol=1e-8);

propertynames(Achol)
#output:  (:U, :L, :p, :P)

Achol.rank
#output: 5

propertynames(Achol, true)
#output: (:U, :L, :p, :P, :factors, :uplo, :piv, :rank, :tol, :info)


getproperty(Achol, :rank)
#output: 5

My question: Why is the rank a private field for the pivoted cholesky?

2

This is the rank return value of the LAPACK pstrf function (Cholesky with complete pivoting).

I’m not sure how reliable a rank estimate this is? LU can be an unreliable way to compute rank, even with complete pivoting (see e.g. the links in this discussion), but I’m not sure if complete-pivoted Cholesky has the same issues.

In any case, it would be good if you could simply call rank(Achol) on the CholeskyPivoted object to access this, similar to how one can call rank on the pivoted QR factorization to get a QR-based rank estimate. Want to file an issue, or better yet a PR?

3

@stevengj Thanks! I think the rank(Achol) function already exists, so an issue (or PR) is not needed for that. I was only trying to understand why this is not a public property.

Did you mean I should submit an issue (or PR) to show rank in the list of symbols when propertynames(Achol) is called? If so I would be happy to give a PR a shot!

4

Oh right, and it is even documented in the CholeskyPivoted docstring. So nothing to do.

I think that since there is a documented API, and it is the standard rank function that is used for other matrices and factorizations, that is preferred to accessing internal fields (which might conceivably change).

5

It’s not the SVD and the error bounds do have factors that grow exponentially with rank, but it’s actually reliable enough in practice that you can consider it to be stable at the same level as QR factorization with column pivoting. (And in fact Cholesky with diagonal pivoting on A^T A corresponds to QR with column pivoting on A). This paper gives an analysis of truncating a pivoted Cholesky factorization. It also includes material on QR factorization with column pivoting and LU with complete pivoting. They are all sort of in the same league of rank determination algorithms that have annoying exponential factors in their error bounds but for which the quantities the exponential factors bound are typically small in practice except on carefully constructed examples. It’s analogous to observing exponential element growth in Gaussian elimination with partial pivoting.

1 Like
6

In general, all fields of objects are internal and subject to change without notice. Unless explicitly documented, that is. But I would never access an internal field of an object, even when documented as public, unless strictly necessary, and even then only grudgingly.

Public fields in Julia always seem out of place and un-idiomatic (imo), since the language nearly universally promotes functional interfaces to objects. If in this case there is a working rank function, it makes perfect sense to hide the internal rank field to discourage field access.

7

Thanks for sharing the reference.