How to quickly do Cholesky factorization for many small matrices?

4 ビュー (過去 30 日間)
Stephan Orzada
Stephan Orzada 2021 年 7 月 23 日
コメント済み: Stephan Orzada 2021 年 8 月 1 日
In my project, I have to check many small matrices for positive-semidefiniteness (PSDness). I use chol() since it is much faster than using eig() to check for negative eigenvalues. Still it is very slow and the most time consuming part of my calculations, since I have to check 10e6 32x32 matrices over and over again. I already tried to use parfor and batch functions but it doesn't work (I also read this in the forums). However, it occured to me that the processor utilization is very low, so I checked what happens when I run the code on two instances of matlab on the same computer. The processor utilization doubled and each script finished in the same time it would have taken if they had been running alone. So it is possible to let chol() run in parallel on the same system. Any ideas how to do this in the same instance of matlab without incurring the problems you get with chol() in a parfor loop?
  12 件のコメント
Bruno Luong
Bruno Luong 2021 年 7 月 31 日
編集済み: Bruno Luong 2021 年 7 月 31 日
In real life 0 eigen value does not exist exactly. CHOL, EIG, EIGS all might fail due to round off.
And Stephan Orzada end up programming his own Cholesky decomposition for his own need, not MATLAB CHOL.
Stephan Orzada
Stephan Orzada 2021 年 8 月 1 日
You are right. For my purpose, I don't care whether it is positive semidefinite or positive definite. Mathematically there is a difference, but it is negligible in many practical cases.

サインインしてコメントする。

回答 (1 件)

Bruno Luong
Bruno Luong 2021 年 7 月 23 日
  2 件のコメント
Stephan Orzada
Stephan Orzada 2021 年 7 月 23 日
編集済み: Stephan Orzada 2021 年 7 月 23 日
Thank you, but unfortunately this doesn't work, since chol in MMX ignores the imaginary part of the complex matrices, and anyway does not support the flag which I need for determining PSDness. I couldn't even tell from the results, because even if the matrix is not PSD (for which the Cholesky decomposition doesn't work) it produces a "result" without any error messages.
Bruno Luong
Bruno Luong 2021 年 7 月 24 日
編集済み: Bruno Luong 2021 年 7 月 24 日
Indeed mmx only deal with real matrix.
But if you are willing to modify the code, you might change the function dpotrf in line 284 of mmx.cpp to zpotrf
Of course you have to take care of retriving MATLAB complex internal interleaved data.
You might ask the author if he can gives you a hand for such task.

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeLinear Algebra についてさらに検索

製品


リリース

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by