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

fixed cholesky error in simulate #314

Merged

Conversation

dietzemarina
Copy link
Contributor

An error occurs when performing cholesky in 'simulate' functions.
When there is at least one element of the diagonal of sys.Q equal to 0.0, the matrix will not be positive definite and the cholesky decomposition won't work.
To bypass this issue, when sys.Q is positive semidefinte (i.e., there are elements of the diagonal of sys.Q equal to 0.0), we add a small number and perform cholesky without errors. Then, we round the elements of the decomposed matrix to recover the original 0.0 values.

@codecov
Copy link

codecov bot commented Jul 29, 2022

Codecov Report

Merging #314 (25f0ca3) into master (407f051) will decrease coverage by 0.23%.
The diff coverage is 40.00%.

❗ Current head 25f0ca3 differs from pull request most recent head 2cfd9d8. Consider uploading reports for the commit 2cfd9d8 to get more accurate results

@@            Coverage Diff             @@
##           master     #314      +/-   ##
==========================================
- Coverage   91.01%   90.77%   -0.24%     
==========================================
  Files          35       35              
  Lines        3296     3308      +12     
==========================================
+ Hits         3000     3003       +3     
- Misses        296      305       +9     
Impacted Files Coverage Δ
src/systems.jl 79.43% <35.71%> (-5.07%) ⬇️
src/models/basicstructural_explanatory.jl 99.00% <100.00%> (ø)

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

Comment on lines 196 to 198
if isposdef(sys.Q[1])
chol_Q = cholesky(sys.Q[1])
elseif ispossemdef(sys.Q[1])
Copy link
Member

Choose a reason for hiding this comment

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

Before this change, if for any reason Q is not positive semidefinite, we would get an error at cholesky which would be clear to the user, but now it will just skip the cholesky step entirely which might result in more confusing errors. Might be good to either change this elseif to else, or have a fallback that throws an error message.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

What do you think about adding an else that gives an error if Q is not positive definite or semidefinite?

Copy link
Member

Choose a reason for hiding this comment

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

Yeah that makes sense to me!

Comment on lines 201 to 203
chol_Q.L[:, :] = round.(chol_Q.L; digits = 6)
chol_Q.U[:, :] = round.(chol_Q.U; digits = 6)
chol_Q.UL[:, :] = round.(chol_Q.UL; digits = 6)
Copy link
Member

Choose a reason for hiding this comment

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

Do we have reasons to be concerned about positive numerical imprecisions that can be magnified in the estimation process because of this rounding? Can we use more digits perhaps?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I'll perform some more tests to check if this brings an important impact. Initially I did this rounding just to be consistent to the original Q matrix, but I guess it might not affect the results significantly.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The process of rounding only affects the estimated scenarios/quantiles. When ingnoring it, there are some differences in the calculated quantiles. What do you think about just increasing the number of digits, as you suggested before? Maybe 10 digits?

Copy link
Member

Choose a reason for hiding this comment

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

Increasing the digits could be enough, let's try 10

@raphaelsaavedra
Copy link
Member

Thanks 🙂

@raphaelsaavedra
Copy link
Member

Can you also bump the package version so we can make a release with this change? Thanks

Copy link
Member

@raphaelsaavedra raphaelsaavedra left a comment

Choose a reason for hiding this comment

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

Thanks!

@raphaelsaavedra raphaelsaavedra merged commit 1dfc0c2 into LAMPSPUC:master Aug 4, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants