Hi, folks. I've run into a number of problems with matrix solution routines in GiNaC, and I'd like to sort them out one by one. The basic overview of the problems is this: 1) matrix::gauss_elimination checks for zero entries via a simple .is_zero() call without a .normal() call before, so it misses zeros in non-normal form. For example it will not detect that this expression: 1/(x - 1) + 1/(x + 1) - 2*x/(x*x - 1) ... is zero. This leads to "division by zero" errors. Test case: #include <ginac/ginac.h> using namespace GiNaC; int main() { symbol x("x"); matrix mx { {1,0,0}, {0,1/(x+1)-1/(x-1)+2/(x*x-1),1}, {0,1,2}, }; matrix zero(mx.rows(), 1); matrix tmp(mx.cols(), 1); for (unsigned i = 0; i < tmp.nops(); i++) { tmp.let_op(i) = symbol(); } mx.solve(tmp, zero, solve_algo::gauss); return 0; } 2) matrix matrix::solve does the same erorr, so if gauss_elimination returned a matrix with un-normal zeros, it too may error out with "division by zero". Test case: #include <ginac/ginac.h> using namespace GiNaC; int main() { symbol x("x"); matrix mx { {1,0,0}, {0,1/(x+1)-1/(x-1)+2/(x*x-1),1}, {0,0,0}, }; matrix zero(mx.rows(), 1); matrix tmp(mx.cols(), 1); for (unsigned i = 0; i < tmp.nops(); i++) { tmp.let_op(i) = symbol(); } mx.solve(tmp, zero, solve_algo::gauss); return 0; } 3) The default echelon form algorithm, Bareiss elimination, (matrix::fraction_free_elimination) is unbearably slow for sparse matrices, and should never be used with them. This is because of explosive common factor growth during its evaluation. For an illustration, using matrix::fraction_free_elimination on a matrix like this: {{a11,a12,a13,a14}, { 0,a22,a23,a24}, { 0, 0,a33,a34}, { 0, 0, 0,a44}} ... will give you a result like this: {{a11, a12, a13, a14}, { 0,a11*a22, a11*a23, a11*a24}, { 0, 0,a11*a22*a33, a11*a22*a34}, { 0, 0, 0,a11*a22*a33*a44}} Now imagine a similar 100x100 matrix with slightly more complex diagonal elements. Now, there are multiple ways of dealing with these problems. The practical workaround for the moment is to: a) use solve_algo::gauss for all sparse-ish matrices b) normal()ize the whole matrix before passing it into solve() There are places where this workaround is not available. One of such places is matrix::inverse(): it doesn't provide an 'algo' argument, and always defaults to solve_algo::automatic (which means fraction_free_elimination). So for the start I propose to add that argument (see the attached patch). The reservation I have with this patch is other functions (solve() and determinant()) have their own enumerations for the 'algo' argument, while inverse() will be reusing the one from solve(). Not the cleanest thing in the world from API usage point of view. What do you think of it?
On Thu, 18 Jan 2018 15:10:39 +0100, Vitaly Magerya <vmagerya@gmail.com> said: VM> 1) matrix::gauss_elimination checks for zero entries via a VM> simple .is_zero() call without a .normal() call before, so it VM> misses zeros in non-normal form.
From time to time I am asking myself either .is_zero() needs to call .normal() first and then make the test. This would look more natural in my opinion. -- Vladimir V. Kisil http://www.maths.leeds.ac.uk/~kisilv/ Book: Geometry of Mobius Transformations http://goo.gl/EaG2Vu Software: Geometry of cycles http://moebinv.sourceforge.net/
On 01/18/2018 03:19 PM, Vladimir V. Kisil wrote:
From time to time I am asking myself either .is_zero() needs to call .normal() first and then make the test. This would look more natural in my opinion.
It will, but usually if you've called normal() on something you don't want to throw the result away, you want to store it. This is because normal() is slow (GCDs and stuff). Auto-normal()ing in is_zero() and throwing away the result is a waste, and can be a major performance hit for some applications. On a related note: when testing for zero you don't really need the full normal() calculation, which involves canceling common divisors of numerator and denominator -- you can keep those in; if you have any in the final result, then the expression is not a zero anyway. Basically, there's a place for normal() equivalent which never calls gcd(). On another related note: there's a flag that tells you if an expression is expanded (status_flag::expanded), which presumably exists so that expand() repeated twice would be a quick no-op. There's no such flag for normal() though, and I'd argue it's much more needed there than in expand(). With such a flag you could pretty much put normal() all over the place without fearing for wasted computation: performance-minded code would have already normal()ized whatever it's passing into is_zero() and the like, so no performance penalty there, and more careless code would get the benefit of always obtaining correct results even when it forgot to normal()ize it's expressions.
On Thu, 18 Jan 2018 15:48:32 +0100, Vitaly Magerya <vmagerya@gmail.com> said:
VM> On 01/18/2018 03:19 PM, Vladimir V. Kisil wrote: >> From time to time I am asking myself either .is_zero() needs to >> call .normal() first and then make the test. This would look more >> natural in my opinion. VM> It will, but usually if you've called normal() on something you VM> don't want to throw the result away, you want to store it. This VM> is because normal() is slow (GCDs and stuff). Auto-normal()ing VM> in is_zero() and throwing away the result is a waste, and can be VM> a major performance hit for some applications. VM> On a related note: when testing for zero you don't really need VM> the full normal() calculation, which involves canceling common VM> divisors of numerator and denominator -- you can keep those in; VM> if you have any in the final result, then the expression is not VM> a zero anyway. Basically, there's a place for normal() VM> equivalent which never calls gcd(). VM> On another related note: there's a flag that tells you if an VM> expression is expanded (status_flag::expanded), which presumably VM> exists so that expand() repeated twice would be a quick VM> no-op. There's no such flag for normal() though, and I'd argue VM> it's much more needed there than in expand(). With such a flag VM> you could pretty much put normal() all over the place without VM> fearing for wasted computation: performance-minded code would VM> have already normal()ized whatever it's passing into is_zero() VM> and the like, so no performance penalty there, and more careless VM> code would get the benefit of always obtaining correct results VM> even when it forgot to normal()ize it's expressions. I agree with all your points. -- Vladimir V. Kisil http://www.maths.leeds.ac.uk/~kisilv/ Book: Geometry of Mobius Transformations http://goo.gl/EaG2Vu Software: Geometry of cycles http://moebinv.sourceforge.net/
Hi Vitaly, On 18.01.2018 15:48, Vitaly Magerya wrote:
It will, but usually if you've called normal() on something you don't want to throw the result away, you want to store it. This is because normal() is slow (GCDs and stuff). Auto-normal()ing in is_zero() and throwing away the result is a waste, and can be a major performance hit for some applications.
On a related note: when testing for zero you don't really need the full normal() calculation, which involves canceling common divisors of numerator and denominator -- you can keep those in; if you have any in the final result, then the expression is not a zero anyway. Basically, there's a place for normal() equivalent which never calls gcd().
It's better to keep .is_zero() lightweight. Adding clever things not only makes it slow, it also might be a slippery road: Where should one stop? How clever should we make it? (Cf. Richardson's theorem.)
On another related note: there's a flag that tells you if an expression is expanded (status_flag::expanded), which presumably exists so that expand() repeated twice would be a quick no-op. There's no such flag for normal() though, and I'd argue it's much more needed there than in expand(). With such a flag you could pretty much put normal() all over the place without fearing for wasted computation: performance-minded code would have already normal()ized whatever it's passing into is_zero() and the like, so no performance penalty there, and more careless code would get the benefit of always obtaining correct results even when it forgot to normal()ize it's expressions.
This sounds very reasonable. Would you like to try your hand at a patch? All my best, -richy.
On 01/25/2018 10:10 PM, Richard B. Kreckel wrote:
On another related note: there's a flag that tells you if an expression is expanded (status_flag::expanded), which presumably exists so that expand() repeated twice would be a quick no-op. There's no such flag for normal() though, and I'd argue it's much more needed there than in expand(). With such a flag you could pretty much put normal() all over the place without fearing for wasted computation: performance-minded code would have already normal()ized whatever it's passing into is_zero() and the like, so no performance penalty there, and more careless code would get the benefit of always obtaining correct results even when it forgot to normal()ize it's expressions.
This sounds very reasonable. Would you like to try your hand at a patch?
I've been testing the attached patch for a bit, but I do not see as big of a performance gain on my workload as I hoped for. This is because normalizing an already normal expression is already much, much faster than normalizing it in the first place, so this optimization is nice, but not as big of a deal. I haven't covered all the places status_flag::normal could be exploited though, so don't commit it just yet.
On 01/18/2018 03:10 PM, I wrote:
[...] Now, there are multiple ways of dealing with these problems. The practical workaround for the moment is to: a) use solve_algo::gauss for all sparse-ish matrices b) normal()ize the whole matrix before passing it into solve()
There are places where this workaround is not available. One of such places is matrix::inverse()
The second place where a better algorithm is needed is matrix::rank(). Currently it blindly fraction_free_elimination(), which should never be used with sparse matrices. BTW, I noticed that matrix::determinant() has some logic to determine matrix sparseness, and it would actually switch away from fraction_free_elimination() for sparse matrices (it will use minor expansion method instead), where sparseness is determined by comparing the number of zero cells to the number of all cells divided by 5. I think it would be beneficial if the same sort of logic would be used in matrix::solve() too. (See the attached patch for my proposal). This would additionally alleviate the need for an 'algo' argument in matrix::inverse(), since solve_algo::automatic would now perform the sensible thing for sparse matrices (which is Gauss elimination). This of course doesn't solve the problem with matrix::rank(). I think there are two good ways to proceed: 1) Apply the attached patch, and make matrix::rank() use similar logic to determine which *_elimination() to call. 2) Apply the attached patch, and add 'algo' arguments to both matrix::inverse() and matrix::rank() (probably using separate enumerations for each).
Hi Vitaly, On 18.01.2018 15:10, Vitaly Magerya wrote:
Hi, folks. I've run into a number of problems with matrix solution routines in GiNaC, and I'd like to sort them out one by one.
The basic overview of the problems is this:
1) matrix::gauss_elimination checks for zero entries via a simple .is_zero() call without a .normal() call before, so it misses zeros in non-normal form. For example it will not detect that this expression:
1/(x - 1) + 1/(x + 1) - 2*x/(x*x - 1)
... is zero. This leads to "division by zero" errors.
Test case:
#include <ginac/ginac.h> using namespace GiNaC; int main() { symbol x("x"); matrix mx { {1,0,0}, {0,1/(x+1)-1/(x-1)+2/(x*x-1),1}, {0,1,2}, }; matrix zero(mx.rows(), 1); matrix tmp(mx.cols(), 1); for (unsigned i = 0; i < tmp.nops(); i++) { tmp.let_op(i) = symbol(); } mx.solve(tmp, zero, solve_algo::gauss); return 0; }
Notice that the comment of matrix::gauss_elimination says that "the algorithm is ok for matrices with numeric coefficients but quite unsuited for symbolic matrices" and there are similar comments about Gauss elimination inside flags.h. Bareiss is the favored But your point seems to be that this comment is incorrect in the first place, since for a certain level of sparseness, Gauss would perform better. See below... I also notice another inconsistency: matrix::fraction_free_elimination is prepared to cope with non-normal zeros, the other algorithms are not. Maybe, that should be harmonized (by making the other algorithms normal their elements first).
2) matrix matrix::solve does the same erorr, so if gauss_elimination returned a matrix with un-normal zeros, it too may error out with "division by zero".
Test case:
#include <ginac/ginac.h> using namespace GiNaC; int main() { symbol x("x"); matrix mx { {1,0,0}, {0,1/(x+1)-1/(x-1)+2/(x*x-1),1}, {0,0,0}, }; matrix zero(mx.rows(), 1); matrix tmp(mx.cols(), 1); for (unsigned i = 0; i < tmp.nops(); i++) { tmp.let_op(i) = symbol(); } mx.solve(tmp, zero, solve_algo::gauss); return 0; }
Same comment as above.
3) The default echelon form algorithm, Bareiss elimination, (matrix::fraction_free_elimination) is unbearably slow for sparse matrices, and should never be used with them. This is because of explosive common factor growth during its evaluation.
For an illustration, using matrix::fraction_free_elimination on a matrix like this:
{{a11,a12,a13,a14}, { 0,a22,a23,a24}, { 0, 0,a33,a34}, { 0, 0, 0,a44}}
... will give you a result like this:
{{a11, a12, a13, a14}, { 0,a11*a22, a11*a23, a11*a24}, { 0, 0,a11*a22*a33, a11*a22*a34}, { 0, 0, 0,a11*a22*a33*a44}}
Now imagine a similar 100x100 matrix with slightly more complex diagonal elements.
Sure, but doing Gauss elemination is gambling with luck: for dense non-numeric matrices, Gauss is much worse! Do you have an idea how to improve the heuristics?
Now, there are multiple ways of dealing with these problems. The practical workaround for the moment is to: a) use solve_algo::gauss for all sparse-ish matrices b) normal()ize the whole matrix before passing it into solve()
There are places where this workaround is not available. One of such places is matrix::inverse(): it doesn't provide an 'algo' argument, and always defaults to solve_algo::automatic (which means fraction_free_elimination). So for the start I propose to add that argument (see the attached patch). The reservation I have with this patch is other functions (solve() and determinant()) have their own enumerations for the 'algo' argument, while inverse() will be reusing the one from solve(). Not the cleanest thing in the world from API usage point of view. What do you think of it?
No objection to re-using solve_algo. All my best, -richy.
On 01/25/2018 11:10 PM, Richard B. Kreckel wrote:
Notice that the comment of matrix::gauss_elimination says that "the algorithm is ok for matrices with numeric coefficients but quite unsuited for symbolic matrices" and there are similar comments about Gauss elimination inside flags.h. Bareiss is the favored
But your point seems to be that this comment is incorrect in the first place, since for a certain level of sparseness, Gauss would perform better. See below...
Right, I don't trust that comment; it runs contrary to my experience.
I also notice another inconsistency: matrix::fraction_free_elimination is prepared to cope with non-normal zeros, the other algorithms are not. Maybe, that should be harmonized (by making the other algorithms normal their elements first).
This is probably because gauss_elimination was not widely used for symbolic matrices, and thus did not receive enough of hardening.
Sure, but doing Gauss elemination is gambling with luck: for dense non-numeric matrices, Gauss is much worse!
Do you have an idea how to improve the heuristics?
Yeah, my proposal is here: https://www.ginac.de/pipermail/ginac-list/2018-January/002174.html BTW, it was my impression that Gauss elimination is not *that* bad for dense matrices. I'm attaching a test case which constructs a dense matrix with some random-ish polynomials and solves it via both Gauss and Bareiss elimination methods. Here are some typical results: matrix size gauss time bareiss time gauss/bareiss 2x2 0.00416815 0.00361022 1.15454 3x3 0.0195643 0.0141638 1.38129 4x4 0.0896266 0.0666258 1.34522 5x5 0.205633 0.185446 1.10885 6x6 0.730309 0.979617 0.745505 7x7 1.29834 3.14298 0.413093 8x8 2.46307 9.70648 0.253755 9x9 4.4289 (segfaults) 10x10 8.39434 (segfaults) ... It's actually worse than I thought it was. I think I need to reconsider my proposal above; it seems that we need to use Gauss elimination for *all* the matrices.
On 01/26/2018 02:24 PM, I wrote:
matrix size gauss time bareiss time gauss/bareiss 2x2 0.00416815 0.00361022 1.15454 3x3 0.0195643 0.0141638 1.38129 4x4 0.0896266 0.0666258 1.34522 5x5 0.205633 0.185446 1.10885 6x6 0.730309 0.979617 0.745505 7x7 1.29834 3.14298 0.413093 8x8 2.46307 9.70648 0.253755 9x9 4.4289 (segfaults) 10x10 8.39434 (segfaults)
These segfaults have this sort of backtrace: #1 0x00007ffff7a07f79 in GiNaC::ex::construct_from_basic(...) #2 0x00007ffff79e19e2 in GiNaC::basic::subs_one_level(...) const #3 0x00007ffff7b3d235 in GiNaC::power::subs(...) const #4 0x00007ffff7a0ba11 in GiNaC::expairseq::subschildren(...) const #5 0x00007ffff7a0bd64 in GiNaC::expairseq::subs(...) const #6 0x00007ffff7a08d4e in GiNaC::ex::subs(...) const #7 0x00007ffff7aeb6c0 in GiNaC::gcd(...) #8 0x00007ffff7aea806 in GiNaC::gcd_pf_pow(...) #9 0x00007ffff7aebbd0 in GiNaC::gcd(...) #10 0x00007ffff7aeac26 in GiNaC::gcd_pf_pow(...) #11 0x00007ffff7aebbd0 in GiNaC::gcd(...) (repeat the last two lines a whole bunch of times) #20147 0x00007ffff7aebbd0 in GiNaC::gcd(...) #20148 0x00007ffff7aeac26 in GiNaC::gcd_pf_pow(...) #20149 0x00007ffff7aea2f8 in GiNaC::gcd_pf_pow(...) #20150 0x00007ffff7aebbd0 in GiNaC::gcd(...) #20151 0x00007ffff7af214a in GiNaC::frac_cancel(...) #20152 0x00007ffff7af2b13 in GiNaC::mul::normal(...) const #20153 0x00007ffff7af2dd0 in GiNaC::add::normal(...) const #20154 0x00007ffff7af2947 in GiNaC::mul::normal(...) const #20155 0x00007ffff7ae613c in GiNaC::ex::normal() const #20156 0x00007ffff7ac73cc in GiNaC::matrix::solve(...) const #20157 0x0000000000404671 in main () at gauss-vs-bareiss.cpp:35 I think this is the same problem Patrick Schulz reported a bit earlier over here: https://www.ginac.de/pipermail/ginac-list/2018-January/002168.html
Hi, On 28.01.2018 19:35, Alexey Sheplyakov wrote:
27.01.2018, 02:42, "Vitaly Magerya" <vmagerya@gmail.com>:
These segfaults have this sort of backtrace:
#1 0x00007ffff7a07f79 in GiNaC::ex::construct_from_basic(...) #2 0x00007ffff79e19e2 in GiNaC::basic::subs_one_level(...) const #3 0x00007ffff7b3d235 in GiNaC::power::subs(...) const #4 0x00007ffff7a0ba11 in GiNaC::expairseq::subschildren(...) const #5 0x00007ffff7a0bd64 in GiNaC::expairseq::subs(...) const #6 0x00007ffff7a08d4e in GiNaC::ex::subs(...) const #7 0x00007ffff7aeb6c0 in GiNaC::gcd(...) #8 0x00007ffff7aea806 in GiNaC::gcd_pf_pow(...) #9 0x00007ffff7aebbd0 in GiNaC::gcd(...) #10 0x00007ffff7aeac26 in GiNaC::gcd_pf_pow(...) #11 0x00007ffff7aebbd0 in GiNaC::gcd(...)
(repeat the last two lines a whole bunch of times)
#20147 0x00007ffff7aebbd0 in GiNaC::gcd(...) #20148 0x00007ffff7aeac26 in GiNaC::gcd_pf_pow(...) #20149 0x00007ffff7aea2f8 in GiNaC::gcd_pf_pow(...) Looks like gcd() goes to an infinite recursion. Could you please make a breakpoint in GiNaC::gcd_pf_pow and print expressions which cause this misbehavior?
Well, it's not an infinite recursion. But computations like gcd(expand(x^1000*(1+x)^2),x^1000) cause GiNaC to recurse so deeply that it eventually overflows the stack. Also, they take way too long. They appear a lot in Vitaly's example. This goofy recursion ought to be shortcut, somehow. -richy.
On 26.01.2018 23:42, Vitaly Magerya wrote:
On 01/26/2018 02:24 PM, I wrote:
matrix size gauss time bareiss time gauss/bareiss 2x2 0.00416815 0.00361022 1.15454 3x3 0.0195643 0.0141638 1.38129 4x4 0.0896266 0.0666258 1.34522 5x5 0.205633 0.185446 1.10885 6x6 0.730309 0.979617 0.745505 7x7 1.29834 3.14298 0.413093 8x8 2.46307 9.70648 0.253755 9x9 4.4289 (segfaults) 10x10 8.39434 (segfaults)
These segfaults have this sort of backtrace:
#1 0x00007ffff7a07f79 in GiNaC::ex::construct_from_basic(...) #2 0x00007ffff79e19e2 in GiNaC::basic::subs_one_level(...) const #3 0x00007ffff7b3d235 in GiNaC::power::subs(...) const #4 0x00007ffff7a0ba11 in GiNaC::expairseq::subschildren(...) const #5 0x00007ffff7a0bd64 in GiNaC::expairseq::subs(...) const #6 0x00007ffff7a08d4e in GiNaC::ex::subs(...) const #7 0x00007ffff7aeb6c0 in GiNaC::gcd(...) #8 0x00007ffff7aea806 in GiNaC::gcd_pf_pow(...) #9 0x00007ffff7aebbd0 in GiNaC::gcd(...) #10 0x00007ffff7aeac26 in GiNaC::gcd_pf_pow(...) #11 0x00007ffff7aebbd0 in GiNaC::gcd(...)
(repeat the last two lines a whole bunch of times)
#20147 0x00007ffff7aebbd0 in GiNaC::gcd(...) #20148 0x00007ffff7aeac26 in GiNaC::gcd_pf_pow(...) #20149 0x00007ffff7aea2f8 in GiNaC::gcd_pf_pow(...) #20150 0x00007ffff7aebbd0 in GiNaC::gcd(...) #20151 0x00007ffff7af214a in GiNaC::frac_cancel(...) #20152 0x00007ffff7af2b13 in GiNaC::mul::normal(...) const #20153 0x00007ffff7af2dd0 in GiNaC::add::normal(...) const #20154 0x00007ffff7af2947 in GiNaC::mul::normal(...) const #20155 0x00007ffff7ae613c in GiNaC::ex::normal() const #20156 0x00007ffff7ac73cc in GiNaC::matrix::solve(...) const #20157 0x0000000000404671 in main () at gauss-vs-bareiss.cpp:35
I think this is the same problem Patrick Schulz reported a bit earlier over here:
https://www.ginac.de/pipermail/ginac-list/2018-January/002168.html
Stack overflow. Increasing the stack size to 1GiB with ulimit -s 1048576 makes the 9x9 case work for me in about 5 minutes. But there must be something goofy going on here that ought to be debugged. -richy.
Hi, On 26.01.2018 23:42, Vitaly Magerya wrote:
I think this is the same problem Patrick Schulz reported a bit earlier over here:
https://www.ginac.de/pipermail/ginac-list/2018-January/002168.html
Could you guys please re-test against current git HEAD? The crashes you reported should be fixed and timings should have improved. Based on this, it may now make more sense to compare algorithms for solving linear equations. Cheers -richy.
On 01/31/2018 12:26 PM, Richard B. Kreckel wrote:
Could you guys please re-test against current git HEAD? The crashes you reported should be fixed and timings should have improved.
The crashes have vanished in my test case. Well done.
Based on this, it may now make more sense to compare algorithms for solving linear equations.
With today's commits, the same toy benchmark gives these timings (with 5 measurements per matrix size): matrix gauss time/bareiss time size average min max 2x2 1.261 1.181 1.313 3x3 1.569 1.385 1.661 4x4 1.549 1.490 1.647 5x5 1.629 1.598 1.692 6x6 1.409 1.191 1.713 7x7 1.014 0.742 1.335 8x8 0.843 0.385 1.542 9x9 0.414 0.172 0.599 10x10 0.184 0.061 0.622 11x11 0.096 0.029 0.210 So, Bareiss performs much better than previously. Returning to the previous topic of switching to Gauss elimination for sparse matrices: my goal here was to demonstrate that on dense matrices Gauss elimination is at least comparable to Bareiss. I think this point should be uncontroversial now.
On 31.01.2018 17:42, Vitaly Magerya wrote:
Returning to the previous topic of switching to Gauss elimination for sparse matrices: my goal here was to demonstrate that on dense matrices Gauss elimination is at least comparable to Bareiss. I think this point should be uncontroversial now.
I intend to release 1.7.3 soon. If you can prepare a patch with commit message and all after that, it will be applied. Let's also pick up your work on the normalized flag after that. -richy.
On 18.01.2018 15:10, Vitaly Magerya wrote:
3) The default echelon form algorithm, Bareiss elimination, (matrix::fraction_free_elimination) is unbearably slow for sparse matrices, and should never be used with them. This is because of explosive common factor growth during its evaluation.
For an illustration, using matrix::fraction_free_elimination on a matrix like this:
{{a11,a12,a13,a14}, { 0,a22,a23,a24}, { 0, 0,a33,a34}, { 0, 0, 0,a44}}
... will give you a result like this:
{{a11, a12, a13, a14}, { 0,a11*a22, a11*a23, a11*a24}, { 0, 0,a11*a22*a33, a11*a22*a34}, { 0, 0, 0,a11*a22*a33*a44}}
Now imagine a similar 100x100 matrix with slightly more complex diagonal elements.
To consider an echelon input matrix is indeed an instructive example. In the first elimination step, Bareiss fraction free elimination produces a 3x3 submatrix where all elements have the common factor a11: {{a11, a12, a13, a14}, { 0,a11*a22,a11*a23,a11*a24}, { 0, 0,a11*a33,a11*a34}, { 0, 0, 0,a11*a44}} This is because the elements in the first column below the pivot element a11 are all zero. It should be possible to use this observation and simply skip this fraction free elimination step! This would avoid the term growth if the elements below the pivot are all zero. Of course, we mustn't divide by a11 in the next step (because it won't divide). And this won't work for determinant computation (because the determinant ends up in the lower right element with the Bareiss elimination). I'm not against tuning the heuristic choice which algorithm to use. But with this modified fraction free elimination, the input matrix would be equal to the output matrix: No term growth for matrices in echelon form and no GCD computation ever. Cheers -richy.
On 01/28/2018 05:58 PM, Richard B. Kreckel wrote:
In the first elimination step, Bareiss fraction free elimination produces a 3x3 submatrix where all elements have the common factor a11:
{{a11, a12, a13, a14}, { 0,a11*a22,a11*a23,a11*a24}, { 0, 0,a11*a33,a11*a34}, { 0, 0, 0,a11*a44}}
This is because the elements in the first column below the pivot element a11 are all zero.
It should be possible to use this observation and simply skip this fraction free elimination step! This would avoid the term growth if the elements below the pivot are all zero.
Of course, we mustn't divide by a11 in the next step (because it won't divide). And this won't work for determinant computation (because the determinant ends up in the lower right element with the Bareiss elimination). This was my thought as well. Turns out that not all of the factors can be eliminated this way. Consider a block-triangular matrix like this one:
{{a11, a12, a13, a14, a15}, {a21, a22, a23, a24, a25}, { 0, a32, a33, a34, a35}, { 0, 0, 0, a44, a45}, { 0, 0, 0, 0, a55}} If we'd just check for zero column under the pivot element, and skip the elimination operation along with subsequent division, we'll get this result: {{a11, a12, a13, a14, a15}, { 0, ..., ..., ..., ...}, { 0, 0, ..., ..., ...}, { 0, 0, 0, (a22*a11-a12*a21)*a44, (a22*a11-a12*a21)*a45}, { 0, 0, 0, 0, (a22*a11-a12*a21)*a55}} As you can see there's still an additional common factor. Now, I did found a way to remove all such factors by tracking them separately, but this costs 2 GCDs per pivot. It's faster than the original algorithm for sparse matrices, but still not as fast as Gauss, at least in my limited testing. See the attachment for the two modified Bareiss elimination schemes: one that only special-cases the zero column case, and the one that eliminates all the factors. Note: there's a paper titled "Fraction Free Gaussian Elimination for Sparse Matrices" [1], but I haven't yet took the time to see what sort of a modification are they using. [1] http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.51.5429
participants (4)
-
Alexey Sheplyakov
-
Richard B. Kreckel
-
Vitaly Magerya
-
Vladimir V. Kisil