FazBrowse GitHub Viewer | Trending |
URL:
| Home
Tools: [Download Repo ZIP]   [Original HTTPS Page]

MAINT: use cholesky_up gufunc for upper Cholesky decomposition by mhvk · Pull Request #25473 · numpy/numpy · GitHub

/ numpy Public

MAINT: use cholesky_up gufunc for upper Cholesky decomposition - #25473

Merged
mhvk merged 1 commit into
numpy:mainfrom
mhvk:cholesky-upper
Jan 3, 2024
Merged

mhvk merged 1 commit into
numpy:mainfrom
mhvk:cholesky-upper

Conversation

mhvk commented Dec 22, 2023

Copy link
Copy Markdown
Contributor

As I was in gufunc-land anyway, thought I might as well do the follow-up of #25388, of using the actual BLAS routine for upper-triangle cholesky decomposition (rather than .mT.conj(), which uses extra memory).

Fixes #25457

Comment on lines +991 to +1000
{
size_t i, j;
for (i = 0; i < n-1; ++i) {
for (j = i+1; j < n; ++j) {
matrix[j] = numeric_limits<typ>::zero;
}
matrix += n;
}
}

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Choose a reason Spam Abuse Off Topic Outdated Duplicate Resolved Low Quality

I've got a question: I see that triu_matrix got renamed to tril_matrix. Was it incorrect?
AFAIU triu_matrix zeroes the upper part of the matrix instead of the lower (assuming row-major-order), or is it column-major implementation?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Choose a reason Spam Abuse Off Topic Outdated Duplicate Resolved Low Quality

This is just what git diff does - really what happened here is that I added triu_matrix - `tril_matrix continues to live below!

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Choose a reason Spam Abuse Off Topic Outdated Duplicate Resolved Low Quality

Here's an implementation of triu_matrix in the current main branch:

template<typename typ>
static inline void
triu_matrix(typ *matrix, size_t n)
{
size_t i, j;
matrix += n;
for (i = 1; i < n; ++i) {
for (j = 0; j < i; ++j) {
matrix[j] = numeric_limits<typ>::zero;
}
matrix += n;
}
}

And here's triu_matrix from your branch:
https://github.com/mhvk/numpy/blob/502bfb27968746f70155130b5fa9579e4d9a3a85/numpy/linalg/umath_linalg.cpp#L987-L1000

It looks like the implementation of triu_matrix changed here, right?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Choose a reason Spam Abuse Off Topic Outdated Duplicate Resolved Low Quality

What I mean: tril_matrix is the one that is new and triu_matrix is the existing one, and I wonder why the existing one changed implementation. I might be missing something.

mtsokol Jan 2, 2024
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Choose a reason Spam Abuse Off Topic Outdated Duplicate Resolved Low Quality

I think this is a typo in the main branch version. In def cholesky there's assert(uplo == 'L'); but later triu_matrix(params.A, params.N); is called.
So I think the change here is correct (assuming row-major order), but the release note should say that triu_matrix implementation changed and tril_matrix was added.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Choose a reason Spam Abuse Off Topic Outdated Duplicate Resolved Low Quality

Ah, yes, sorry! I did change the implementation, as indeed for a triangular-upper matrix, one should be zeroing the lower part. I now recall puzzling over this too, and changing the comment above to try to make clear what was actually done.

mtsokol Jan 2, 2024
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Choose a reason Spam Abuse Off Topic Outdated Duplicate Resolved Low Quality

Ok! I've got one more question: Here triu_matrix is zeroing upper triangular if one assumes the input is row-major order and lower triangular if one assumes the input column-major order, right?

Then maybe we could have just a docstring that explains why:

triu_matrix(typ *matrix, size_t n)
{
    size_t i, j;
    for (i = 0; i < n-1; ++i) {
        for (j = i+1; j < n; ++j) {
            matrix[j] = numeric_limits<typ>::zero;
        }
        matrix += n;
    }
}

looks like it zeroes upper triangular but the function's name says triu? Just a comment that matrix is column major (I think so).
WDYT?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Choose a reason Spam Abuse Off Topic Outdated Duplicate Resolved Low Quality

Yes, makes sense!

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Choose a reason Spam Abuse Off Topic Outdated Duplicate Resolved Low Quality

OK, added the extra comment (really good to have indeed!)

mtsokol left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Choose a reason Spam Abuse Off Topic Outdated Duplicate Resolved Low Quality

Apart from the one posted question, LGTM. (I think a release note can be added that says that a new ufunc was added)

mhvk commented Jan 2, 2024

Copy link
Copy Markdown
Contributor Author

Do we need a release note even if the ufunc is only used privately? (same for tril_matrix, which is only used inside other code).

mtsokol commented Jan 2, 2024

Copy link
Copy Markdown
Member

Do we need a release note even if the ufunc is only used privately? (same for tril_matrix, which is only used inside other code).

Ah right, I agree there's no need for a note.

mhvk merged commit d906b52 into numpy:main Jan 3, 2024
mhvk deleted the cholesky-upper branch January 3, 2024 01:54
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

ENH: Add a dedicated cholesky_up ufunc

2 participants


Back | FazBrowse Home | New Git URL