-
Notifications
You must be signed in to change notification settings - Fork 56
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
Refactor functions in matrixtools.py #442
base: develop
Are you sure you want to change the base?
Conversation
if not is_hermitian(mx, tol): | ||
return False | ||
if attempt_cholesky: | ||
try: | ||
_ = _spl.cholesky(mx) | ||
return True # Cholesky succeeded | ||
except _spl.LinAlgError: | ||
pass # we fall back on eigenvalue decomposition | ||
evals = _np.linalg.eigvalsh(mx) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The property of being positive semidefinite is generally ascribed to Hermitian matrices, since the value of the quadratic form
The new implementation starts with a check that the matrix is Hermitian. If we do have a Hermitian matrix then the fastest way to check if it's positive definite is a Cholesky decomposition. However, Cholesky will fail on matrices that are only positive semidefinite, so there is no point in trying Cholesky if we expect is_pos_def
will regularly handle singular matrices. I've set the default behavior for attempt_cholesky
to False. If we skip Cholesky or if Cholesky fails then we use Hermitian eigendecomposition.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
A couple of things:
- Should we ensure the docstring is consistent? It says PD in the first line, but return says PSD.
- Both implementations are clearly for PSD. Should the function name be
is_pos_semidef
instead? Pro is that it is more accurate, con is that it is an API change. While this is in matrixtools and therefore could have moderate/high user exposure, it seems like the only use is inis_valid_density_mx
below (which also seems like it has low usage), so maybe the API change would go mostly unnoticed?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In the context of quantum information I find people are lax about positive definite versus positive semidefinite, so I think the function name is fine. But I'll definitely change the docstring to be consistent about really checking that a matrix is PSD.
def remove_dependent_cols(mx, tol=1e-7): | ||
""" | ||
Removes the linearly dependent columns of a matrix. | ||
|
||
Parameters | ||
---------- | ||
mx : numpy.ndarray | ||
The input matrix | ||
|
||
Returns | ||
------- | ||
A linearly independent subset of the columns of `mx`. | ||
""" | ||
last_rank = 0; cols_to_remove = [] | ||
for j in range(mx.shape[1]): | ||
rnk = _np.linalg.matrix_rank(mx[:, 0:j + 1], tol) | ||
if rnk == last_rank: | ||
cols_to_remove.append(j) | ||
else: | ||
last_rank = rnk | ||
#print("Removing %d cols" % len(cols_to_remove)) | ||
return _np.delete(mx, cols_to_remove, axis=1) | ||
|
||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This was called in one function (below), which I reimplemented to rely on independent_columns
.
Get the latest fixes for the ComplementPOVMEffect.to_vector() calls.
Forgot to pull latest develop before attempting the previous merge.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for the great effort @rileyjmurray, as always! Also, I greatly appreciate your preemptive comments on some of the more technical changes, it makes it much easier to follow. :)
There is only one real action item here, for an additional unit test for the attempt_cholesky
codepath in is_pos_def
(and a small docstring update). However, I've also made some non-blocking comments/questions. I 👍'd any of your comments that I read and agreed with but had no active comments or action items for.
# I wanted to remove the implementations of vec and unvec and just in-line equivalent | ||
# code in the few places they were used. However, the fact that they're included in this | ||
# __init__.py file suggests that they might be used outside of pyGSTi itself. | ||
# |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We can check with @kevincyoung and @bruzic3, but I think I'm OK with inlining these in our source, deprecating the functions in pyGSTi, and then assisting/enforcing the relevant changes in the downstream code so that they can be removed if they are inefficient/not best practice.
if not is_hermitian(mx, tol): | ||
return False | ||
if attempt_cholesky: | ||
try: | ||
_ = _spl.cholesky(mx) | ||
return True # Cholesky succeeded | ||
except _spl.LinAlgError: | ||
pass # we fall back on eigenvalue decomposition | ||
evals = _np.linalg.eigvalsh(mx) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
A couple of things:
- Should we ensure the docstring is consistent? It says PD in the first line, but return says PSD.
- Both implementations are clearly for PSD. Should the function name be
is_pos_semidef
instead? Pro is that it is more accurate, con is that it is an API change. While this is in matrixtools and therefore could have moderate/high user exposure, it seems like the only use is inis_valid_density_mx
below (which also seems like it has low usage), so maybe the API change would go mostly unnoticed?
@@ -1372,6 +1424,9 @@ def _findx(a, inds, always_copy=False): | |||
return a_inds | |||
|
|||
|
|||
# TODO: reevaluate the need for this function. It seems like we could just in-line @ |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Generally in favor of this. Not a requirement for approval for this PR though, can punt to future.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Cool. I'll try to make this change then.
This PR resolves #429. I've provided detailed comments explaining my changes to
matrixtools.py
. I've also removed some commented-out code.