Eigen solver for sparse matrix product

I want to use Eigen to compute L_1^-1L_2, where both L_1 and L_2 are lower triangular matrices and are stored as column oriented sparse matrices in Eigen.
I tried the Eigen triangular solver. However, that requires L_2 to be dense.

asked Nov 13 '18 at 22:33

Ting

286

The inverse of a sparse matrix is usually dense, therefore L_1^-1*L_2 will be dense as well. Are you sure, you need to store the matrix, instead of just computing products/solving systems on the fly?

– chtz
Nov 14 '18 at 8:11

@chtz Thanks for your response. I actually have guarantee that L_1^-1 is also sparse.

– Ting
Nov 14 '18 at 16:42

add a comment |

asked Nov 13 '18 at 22:33

Ting

286

The inverse of a sparse matrix is usually dense, therefore L_1^-1*L_2 will be dense as well. Are you sure, you need to store the matrix, instead of just computing products/solving systems on the fly?

– chtz
Nov 14 '18 at 8:11

@chtz Thanks for your response. I actually have guarantee that L_1^-1 is also sparse.

– Ting
Nov 14 '18 at 16:42

add a comment |

asked Nov 13 '18 at 22:33

Ting

286

c++ linear-algebra eigen3

asked Nov 13 '18 at 22:33

Ting

286

asked Nov 13 '18 at 22:33

Ting

286

asked Nov 13 '18 at 22:33

Ting

286

asked Nov 13 '18 at 22:33

Ting

286

asked Nov 13 '18 at 22:33

Ting

286

The inverse of a sparse matrix is usually dense, therefore L_1^-1*L_2 will be dense as well. Are you sure, you need to store the matrix, instead of just computing products/solving systems on the fly?

– chtz
Nov 14 '18 at 8:11

@chtz Thanks for your response. I actually have guarantee that L_1^-1 is also sparse.

– Ting
Nov 14 '18 at 16:42

add a comment |

The inverse of a sparse matrix is usually dense, therefore L_1^-1*L_2 will be dense as well. Are you sure, you need to store the matrix, instead of just computing products/solving systems on the fly?

– chtz
Nov 14 '18 at 8:11

@chtz Thanks for your response. I actually have guarantee that L_1^-1 is also sparse.

– Ting
Nov 14 '18 at 16:42

The inverse of a sparse matrix is usually dense, therefore L_1^-1*L_2 will be dense as well. Are you sure, you need to store the matrix, instead of just computing products/solving systems on the fly?

– chtz
Nov 14 '18 at 8:11

@chtz Thanks for your response. I actually have guarantee that L_1^-1 is also sparse.

– Ting
Nov 14 '18 at 16:42

add a comment |

1 Answer
1

active

oldest

votes

The solve method in fact is not overloaded for sparse rhs, however you can use the solveInPlace method like so (I did not actually try this):

Eigen::SparseMatrix<double> foo(Eigen::SparseMatrix<double> const& L1, Eigen::SparseMatrix<double> const& L2)

 Eigen::SparseMatrix<double> res = L2;
 L1.triangularView<Eigen::Lower>().solveInPlace(res);
 return res;

Still you should consider if you actually need the full matrix.

answered Nov 15 '18 at 11:25

chtz

6,66511232

Thanks for your response. Indeed I need to compute the inverse of a spd matrix K = LL^T, where L is very sparse. My solution is to call triangular solve twice: L_inv = L.triangularView<Lower>().solve(Identity) and L.transpose().triangularView<Upper>.solve(L_inv). Is this the right way to do it in Eigen?

– Ting
Nov 16 '18 at 14:39

Instead of the second solve, you can compute L_inv.transpose()*L_inv

– chtz
Nov 16 '18 at 20:58

But if L is sparse, isn't L.transpose().triangularView<Upper>.solve(L_inv) faster than L_inv.transpose()*L_inv?

– Ting
Nov 19 '18 at 14:06

I can't think of good examples where a sparse solve would be faster than a sparse product. To find out if this happens in your case you need to benchmark!

– chtz
Nov 19 '18 at 15:17

Thanks for your help!

– Ting
Nov 19 '18 at 21:42

add a comment |

Your Answer

StackExchange.ifUsing("editor", function ()
StackExchange.using("externalEditor", function ()
StackExchange.using("snippets", function ()
StackExchange.snippets.init();
);
);
, "code-snippets");

StackExchange.ready(function()
var channelOptions =
tags: "".split(" "),
id: "1"
;
initTagRenderer("".split(" "), "".split(" "), channelOptions);

StackExchange.using("externalEditor", function()
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled)
StackExchange.using("snippets", function()
createEditor();
);

else
createEditor();

);

function createEditor()
StackExchange.prepareEditor(
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: true,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: 10,
bindNavPrevention: true,
postfix: "",
imageUploader:
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
,
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
);

);

draft saved

draft discarded

StackExchange.ready(
function ()
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f53290521%2feigen-solver-for-sparse-matrix-product%23new-answer', 'question_page');

);

Post as a guest

Name

Required, but never shown

1 Answer
1

active

oldest

votes

1 Answer
1

active

oldest

votes

The solve method in fact is not overloaded for sparse rhs, however you can use the solveInPlace method like so (I did not actually try this):

Eigen::SparseMatrix<double> foo(Eigen::SparseMatrix<double> const& L1, Eigen::SparseMatrix<double> const& L2)

 Eigen::SparseMatrix<double> res = L2;
 L1.triangularView<Eigen::Lower>().solveInPlace(res);
 return res;

Still you should consider if you actually need the full matrix.

answered Nov 15 '18 at 11:25

chtz

6,66511232

Thanks for your response. Indeed I need to compute the inverse of a spd matrix K = LL^T, where L is very sparse. My solution is to call triangular solve twice: L_inv = L.triangularView<Lower>().solve(Identity) and L.transpose().triangularView<Upper>.solve(L_inv). Is this the right way to do it in Eigen?

– Ting
Nov 16 '18 at 14:39

Instead of the second solve, you can compute L_inv.transpose()*L_inv

– chtz
Nov 16 '18 at 20:58

But if L is sparse, isn't L.transpose().triangularView<Upper>.solve(L_inv) faster than L_inv.transpose()*L_inv?

– Ting
Nov 19 '18 at 14:06

I can't think of good examples where a sparse solve would be faster than a sparse product. To find out if this happens in your case you need to benchmark!

– chtz
Nov 19 '18 at 15:17

Thanks for your help!

– Ting
Nov 19 '18 at 21:42

add a comment |

The solve method in fact is not overloaded for sparse rhs, however you can use the solveInPlace method like so (I did not actually try this):

Eigen::SparseMatrix<double> foo(Eigen::SparseMatrix<double> const& L1, Eigen::SparseMatrix<double> const& L2)

 Eigen::SparseMatrix<double> res = L2;
 L1.triangularView<Eigen::Lower>().solveInPlace(res);
 return res;

Still you should consider if you actually need the full matrix.

answered Nov 15 '18 at 11:25

chtz

6,66511232

Thanks for your response. Indeed I need to compute the inverse of a spd matrix K = LL^T, where L is very sparse. My solution is to call triangular solve twice: L_inv = L.triangularView<Lower>().solve(Identity) and L.transpose().triangularView<Upper>.solve(L_inv). Is this the right way to do it in Eigen?

– Ting
Nov 16 '18 at 14:39

Instead of the second solve, you can compute L_inv.transpose()*L_inv

– chtz
Nov 16 '18 at 20:58

But if L is sparse, isn't L.transpose().triangularView<Upper>.solve(L_inv) faster than L_inv.transpose()*L_inv?

– Ting
Nov 19 '18 at 14:06

I can't think of good examples where a sparse solve would be faster than a sparse product. To find out if this happens in your case you need to benchmark!

– chtz
Nov 19 '18 at 15:17

Thanks for your help!

– Ting
Nov 19 '18 at 21:42

add a comment |

The solve method in fact is not overloaded for sparse rhs, however you can use the solveInPlace method like so (I did not actually try this):

Eigen::SparseMatrix<double> foo(Eigen::SparseMatrix<double> const& L1, Eigen::SparseMatrix<double> const& L2)

 Eigen::SparseMatrix<double> res = L2;
 L1.triangularView<Eigen::Lower>().solveInPlace(res);
 return res;

Still you should consider if you actually need the full matrix.

answered Nov 15 '18 at 11:25

chtz

6,66511232

The solve method in fact is not overloaded for sparse rhs, however you can use the solveInPlace method like so (I did not actually try this):

Eigen::SparseMatrix<double> foo(Eigen::SparseMatrix<double> const& L1, Eigen::SparseMatrix<double> const& L2)

 Eigen::SparseMatrix<double> res = L2;
 L1.triangularView<Eigen::Lower>().solveInPlace(res);
 return res;

Still you should consider if you actually need the full matrix.

answered Nov 15 '18 at 11:25

chtz

6,66511232

answered Nov 15 '18 at 11:25

chtz

6,66511232

answered Nov 15 '18 at 11:25

chtz

6,66511232

answered Nov 15 '18 at 11:25

chtz

6,66511232

Thanks for your response. Indeed I need to compute the inverse of a spd matrix K = LL^T, where L is very sparse. My solution is to call triangular solve twice: L_inv = L.triangularView<Lower>().solve(Identity) and L.transpose().triangularView<Upper>.solve(L_inv). Is this the right way to do it in Eigen?

– Ting
Nov 16 '18 at 14:39

Instead of the second solve, you can compute L_inv.transpose()*L_inv

– chtz
Nov 16 '18 at 20:58

But if L is sparse, isn't L.transpose().triangularView<Upper>.solve(L_inv) faster than L_inv.transpose()*L_inv?

– Ting
Nov 19 '18 at 14:06

I can't think of good examples where a sparse solve would be faster than a sparse product. To find out if this happens in your case you need to benchmark!

– chtz
Nov 19 '18 at 15:17

Thanks for your help!

– Ting
Nov 19 '18 at 21:42

add a comment |

Thanks for your response. Indeed I need to compute the inverse of a spd matrix K = LL^T, where L is very sparse. My solution is to call triangular solve twice: L_inv = L.triangularView<Lower>().solve(Identity) and L.transpose().triangularView<Upper>.solve(L_inv). Is this the right way to do it in Eigen?

– Ting
Nov 16 '18 at 14:39

Instead of the second solve, you can compute L_inv.transpose()*L_inv

– chtz
Nov 16 '18 at 20:58

But if L is sparse, isn't L.transpose().triangularView<Upper>.solve(L_inv) faster than L_inv.transpose()*L_inv?

– Ting
Nov 19 '18 at 14:06

I can't think of good examples where a sparse solve would be faster than a sparse product. To find out if this happens in your case you need to benchmark!

– chtz
Nov 19 '18 at 15:17

Thanks for your help!

– Ting
Nov 19 '18 at 21:42

Thanks for your response. Indeed I need to compute the inverse of a spd matrix K = LL^T, where L is very sparse. My solution is to call triangular solve twice: L_inv = L.triangularView<Lower>().solve(Identity) and L.transpose().triangularView<Upper>.solve(L_inv). Is this the right way to do it in Eigen?

– Ting
Nov 16 '18 at 14:39

Instead of the second solve, you can compute L_inv.transpose()*L_inv

– chtz
Nov 16 '18 at 20:58

But if L is sparse, isn't L.transpose().triangularView<Upper>.solve(L_inv) faster than L_inv.transpose()*L_inv?

– Ting
Nov 19 '18 at 14:06

I can't think of good examples where a sparse solve would be faster than a sparse product. To find out if this happens in your case you need to benchmark!

– chtz
Nov 19 '18 at 15:17

Thanks for your help!

– Ting
Nov 19 '18 at 21:42

add a comment |

draft saved

draft discarded

Thanks for contributing an answer to Stack Overflow!

Please be sure to answer the question. Provide details and share your research!

But avoid …

Asking for help, clarification, or responding to other answers.

Making statements based on opinion; back them up with references or personal experience.

To learn more, see our tips on writing great answers.

draft saved

draft discarded

Post as a guest

Name

Required, but never shown

Name

Required, but never shown

Name

Required, but never shown

This page is only for reference, If you need detailed information, please check here

搜尋此網誌

Odtnhj