Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Deprecate chol(A, Val{:U/:L}) and move documentation to method definitions #13680

Merged
merged 1 commit into from
Oct 22, 2015

Conversation

andreasnoack
Copy link
Member

It was resently brought up on the mailing list that the syntax for getting the lower triangular version of the Cholesky, chol(A, Val{:L}), is "ugly and confusing".

In 0.3 the syntax was chol(A,:L) which is simpler, but because the Triangular type was split into LowerTriangular and UpperTriangular, it was necessary to introduce two different method signatures for creating the Lower and the Upper version of the Cholesky factor. This let to JuliaLang/LinearAlgebra.jl#268.

In this PR I follow the idea from 3. in JuliaLang/LinearAlgebra.jl#268 and remove the second argument from chol. The reason for this is that the Lower and the Upper versions are in fact the exact same factorization, A=◣*◥, and the only difference is if the lower or the upper part is used for storage. Therefore, it is only a matter of transposing the result to get the other version and if this is done as part of a multiplication of solve, then this is essensially free.

@tkelman tkelman added the linear algebra Linear algebra label Oct 20, 2015
@tkelman
Copy link
Contributor

tkelman commented Oct 20, 2015

-1 on this. Defaulting to upper storage always felt like the wrong default in matlab. I could maybe get behind switching from Val{:U} to UpperTriangular as a flag, but you're bringing in an unrelated type to this API just to use it as a flag. Definitely don't swap the argument order of chol!, that makes no sense since you aren't mutating the type, you're just using it as an option flag.

@mschauer
Copy link
Contributor

Once there are free transposes this would not seem to be a big problem and the physicists and numerical mathematicians convention of conjugating (and transposing) the first factor agrees with the fact that we as well only have one - the physicists' - dot product. On the other hand especially in stochastic applications the lower triangular factor is the more natural object, where variance is defined as E X X'. Hard to call, but in the end, chol(A, Val{:L}) can somehow only be confusing to those who need the lower factor, right?

@andreasnoack
Copy link
Member Author

Defaulting to upper storage always felt like the wrong default in matlab.

Maybe, but it doesn't matter which one you choose and therefore it is not worth the trouble to change the default. It also happens that MATLAB and R use the upper version. For the fun of it, I did some archeology, and it is a very old decision in Julia. See bb1ae75#diff-af7f0b18d31ef100c07d8bf23b112ff6R115. It's from the time when files with Julia code were in the root folder and had extension j.

I think it is better to simplify the method instead of giving unnecessary options.

The rest of your comment is about unexported methods. For unexported methods, the balance between cleaning up and avoiding breakages should be quite different. Using UpperTriangular instead of Val{:U} is more clear and when using the return type for dispatch, I think it makes sense to mimic convert and have the return type argument as the first argument.

@tkelman
Copy link
Contributor

tkelman commented Oct 20, 2015

when using the return type for dispatch, I think it makes sense to mimic convert and have the return type argument as the first argument.

Important distinction - convert isn't mutating. Despite chol! being unexported, it is important in a number of applications. We've been adjusting this internal API repeatedly - swapping the argument order on top of changing how the flags work is even more churn for negligible benefit. In my eyes chol!(UpperTriangular, A) is objectively worse than chol!(A, UpperTriangular), ignoring the existing API.

Or comparing to the other closely related API's around chol!, like cholfact! and the Cholesky constructor, all of them take matrix first, uplo flag second.

@dmbates
Copy link
Member

dmbates commented Oct 20, 2015

I think the reason for using the upper triangular factor in statistical computing in general and in S, the predecessor language to R, is because the triangular matrix R in a QR decomposition of a tall, thin model matrix X is then the same as the Cholesky factor of X'X. Statisticians are taught to think in terms of X, X'X, R'R, etc.

@andreasnoack andreasnoack force-pushed the anj/chol branch 2 times, most recently from b0a26cc to 8a0a6bf Compare October 20, 2015 18:37
@tkelman
Copy link
Contributor

tkelman commented Oct 22, 2015

I can live with this version.

tkelman added a commit that referenced this pull request Oct 22, 2015
Deprecate chol(A, Val{:U/:L}) and move documentation to method definitions
@tkelman tkelman merged commit e102efa into master Oct 22, 2015
@tkelman tkelman deleted the anj/chol branch October 22, 2015 01:26
if uplo == :U
Cholesky(chol!(A, UpperTriangular).data, 'U')
else
Cholesky(chol!(A, UpperTriangular).data, 'U')
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Whoops, was this supposed to be Cholesky(chol!(A, LowerTriangular).data, 'L') ? Definitely missing test coverage if so.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch. I'm really surprised that it is not covered. I'll fix that now.

@andreasnoack
Copy link
Member Author

I've just learned that LINPACK only supports the upper triangular version of the Cholesky which could be the reason why MATLAB and R decided to use the U'U version.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants