Add LDLT Cholesky decomposition #4130
There are no files selected for viewing
This file contains 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
This file contains 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
This file contains 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
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -106,6 +106,45 @@ namespace shogun | |
| return sg_linalg->get_cpu_backend(); | ||
| } | ||
|
|
||
| /** Infer the appropriate backend for linalg operations | ||
| * from the input SGVector or SGMatrix (Container). | ||
| * Raise error if the backends of the Containers conflict. | ||
| * | ||
| * @param a The first SGVector/SGMatrix | ||
| * @param b The second SGVector/SGMatrix | ||
| * @param c The third SGVector/SGMatrix | ||
| * @return @see LinalgBackendBase pointer | ||
| */ | ||
| template <typename T, template <typename> class Container> | ||
| LinalgBackendBase* infer_backend( | ||
| const Container<T>& a, const Container<T>& b, const Container<T>& c) | ||
| { | ||
| if (a.on_gpu() && b.on_gpu() && c.on_gpu()) | ||
| { | ||
| if (sg_linalg->get_gpu_backend()) | ||
| return sg_linalg->get_gpu_backend(); | ||
| else | ||
| { | ||
| SG_SERROR( | ||
| "Vector or matrix is on GPU but no GPU backend registered. \ | ||
| This can happen if the GPU backend was de-activated \ | ||
| after memory has been transferred to GPU.\n"); | ||
| return NULL; | ||
| } | ||
| } | ||
| else if (a.on_gpu() || b.on_gpu() || c.on_gpu()) | ||
| { | ||
| SG_SERROR( | ||
| "Cannot operate with first vector/matrix on_gpu flag(%d),\ | ||
| second vector/matrix on_gpu flag (%d) and third vector/matrix \ | ||
| on_gpu flag (%d).\n", | ||
| a.on_gpu(), b.on_gpu(), c.on_gpu()); | ||
| return NULL; | ||
| } | ||
| else | ||
| return sg_linalg->get_cpu_backend(); | ||
| } | ||
|
|
||
| /** | ||
| * Transfers data to GPU memory. | ||
| * Shallow-copy of SGVector with vector on CPU if GPU backend not | ||
|
|
@@ -590,6 +629,71 @@ namespace shogun | |
| ->cholesky_solver(L, b, lower); | ||
| } | ||
|
|
||
| /** | ||
| * Compute the LDLT cholesky decomposition \f$A = P^{T} L D L^{*} P\f$ | ||
| * or \f$A = P^{T} U^{*} D U P\f$ | ||
| * of a positive semidefinite or negative semidefinite Hermitan matrix | ||
| * | ||
| * @param A The matrix whose LDLT cholesky decomposition is to be | ||
| * computed | ||
| * @param L The matrix that saves the triangular LDLT | ||
| * Cholesky factorization (default: lower) | ||
| * @param d The vector that saves the diagonal of the diagonal matrix D | ||
| * @param p The vector that saves the permutation matrix P as a | ||
| * transposition sequence | ||
| * @param lower Whether to use L as the upper or lower triangular | ||
| * Cholesky factorization (default:lower) | ||
| */ | ||
| template <typename T> | ||
| void ldlt_factor( | ||
| const SGMatrix<T>& A, SGMatrix<T>& L, SGVector<T>& d, | ||
| SGVector<index_t>& p, const bool lower = true) | ||
| { | ||
| REQUIRE( | ||
|
There was a problem hiding this comment. @lisitsyn we need some systematic and compact way to do these checks. The amount of lines of code for such error messages is crazy |
||
| A.num_rows == A.num_cols, | ||
| "Matrix dimensions (%dx%d) are not square\n", A.num_rows, | ||
| A.num_cols); | ||
| REQUIRE( | ||
| A.num_rows == L.num_rows && A.num_cols == L.num_cols, | ||
| "Shape of matrix A (%d, %d) doesn't match matrix L (%d, %d)\n", | ||
| A.num_rows, A.num_cols, L.num_rows, L.num_rows); | ||
| REQUIRE( | ||
| d.vlen == A.num_rows, "Length of vector d (%d) doesn't match " | ||
| "length of diagonal of matrix L (%d)\n", | ||
| d.vlen, A.num_rows); | ||
| REQUIRE( | ||
| p.vlen = A.num_rows, "Length of transpositions vector p (%d) " | ||
| "doesn't match length of diagonal of " | ||
| "matrix L (%d)\n", | ||
| p.vlen, A.num_rows); | ||
|
|
||
| infer_backend(A)->ldlt_factor(A, L, d, p, lower); | ||
| } | ||
|
|
||
| /** | ||
| * Solve the linear equations \f$Ax=b\f$, given the LDLT Cholesky | ||
| * factorization of A, | ||
| * where \f$A\f$ is a positive semidefinite or negative semidefinite | ||
| * Hermitan matrix | ||
|
There was a problem hiding this comment. We could add a |
||
| * | ||
| * @param L Triangular matrix, LDLT Cholesky factorization of A | ||
| * @param d The diagonal of the diagonal matrix D | ||
| * @param p The permuattion matrix P as a | ||
| * transposition sequence | ||
| * @param b Right-hand side array | ||
| * @param lower Whether to use L as the upper or lower triangular | ||
| * Cholesky factorization (default:lower) | ||
| * @return \f$\x\f$ | ||
| */ | ||
| template <typename T> | ||
| SGVector<T> ldlt_solver( | ||
| const SGMatrix<T>& L, const SGVector<T>& d, SGVector<index_t>& p, | ||
| const SGVector<T>& b, const bool lower = true) | ||
| { | ||
| return infer_backend(L, SGMatrix<T>(d), SGMatrix<T>(b)) | ||
| ->ldlt_solver(L, d, p, b, lower); | ||
| } | ||
|
|
||
| /** | ||
| * Vector dot-product that works with generic vectors. | ||
| * | ||
|
|
||
This file contains 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
This file contains 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
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
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.
Can't we make this a bit nicer here, so that any number of things that have the
.on_gpu(), .on_cpu()method implemented is checked?@lisitsyn can you do some magic ?