[GiNaC-list] matrix::solve()-related problems
Vitaly Magerya
vmagerya at gmail.com
Sun Jan 28 19:43:28 CET 2018
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
-------------- next part --------------
// compile and run with:
// c++ bareiss.cpp -lginac -lcln -ldl && ./a.out
#include <ginac/ginac.h>
#include <assert.h>
using namespace GiNaC;
using namespace std;
/* Since the internals of matrix class are protected, we need
* to employ this hack to access them.
*/
class matrix_hack : public matrix {
public:
exvector &get_m() { return m; };
};
void
echelon_form_bareiss_original(matrix &m)
{
unsigned nr = m.rows();
unsigned nc = m.cols();
exvector &mv = ((matrix_hack*)&m)->get_m();
ex z = 1;
unsigned r0 = 0;
for (unsigned c0 = 0; (c0 < nc) && (r0 < nr - 1); c0++) {
unsigned pivot = r0;
// No normalization before is_zero() here, because
// we maintain the matrix normalized throughout the
// algorithm.
while ((pivot < nr) && mv[pivot*nc + c0].is_zero()) pivot++;
if (pivot == nr) {
// The whole column below r0:c0 is zero, let's skip
// it.
continue;
}
if (pivot > r0) {
// Found a non-zero row somewhere below r0; let's
// swap it in.
for (unsigned c = c0; c < nc; c++) {
mv[pivot*nc + c].swap(mv[r0*nc + c]);
}
}
ex a = mv[r0*nc + c0];
for (unsigned r = r0 + 1; r < nr; r++) {
ex b = mv[r*nc + c0];
for (unsigned c = c0 + 1; c < nc; c++) {
ex v;
bool divok = divide(a*mv[r*nc + c] - b*mv[r0*nc + c], z, v);
assert(divok);
// Since v is a polynomial, the call to normal()
// here does not involve GCD, but rather just
// expands sums while preserving common factors,
// which is nomally preferable to full expand().
m[r*nc + c] = normal(v);
}
mv[r*nc + c0] = 0;
}
z = a;
r0++;
}
// Zero out the remaining rows (just in case).
for (unsigned r = r0 + 1; r < nr; r++) {
for (unsigned c = 0; c < nc; c++) {
mv[r*nc + c] = 0;
}
}
}
void
echelon_form_bareiss_skip_zrows(matrix &m)
{
unsigned nr = m.rows();
unsigned nc = m.cols();
exvector &mv = ((matrix_hack*)&m)->get_m();
ex z = 1;
unsigned r0 = 0;
for (unsigned c0 = 0; (c0 < nc) && (r0 < nr - 1); c0++) {
unsigned pivot = r0;
// No normalization before is_zero() here, because
// we maintain the matrix normalized throughout the
// algorithm.
while ((pivot < nr) && mv[pivot*nc + c0].is_zero()) pivot++;
if (pivot == nr) {
// The whole column below r0:c0 is zero, let's skip
// it.
continue;
}
if (pivot > r0) {
// Found a non-zero row somewhere below r0; let's
// swap it in.
for (unsigned c = c0; c < nc; c++) {
mv[pivot*nc + c].swap(mv[r0*nc + c]);
}
}
bool zerocol = true;
for (unsigned r = r0 + 1; r < nr; r++) {
ex b = mv[r*nc + c0];
zerocol = zerocol && b.is_zero();
}
if (zerocol) {
z = 1;
r0++;
continue;
}
ex a = mv[r0*nc + c0];
for (unsigned r = r0 + 1; r < nr; r++) {
ex b = mv[r*nc + c0];
for (unsigned c = c0 + 1; c < nc; c++) {
ex v;
bool divok = divide(a*mv[r*nc + c] - b*mv[r0*nc + c], z, v);
assert(divok);
// Since v is a polynomial, the call to normal()
// here does not involve GCD, but rather just
// expands sums while preserving common factors,
// which is nomally preferable to full expand().
m[r*nc + c] = normal(v);
}
mv[r*nc + c0] = 0;
}
z = a;
r0++;
}
// Zero out the remaining rows (just in case).
for (unsigned r = r0 + 1; r < nr; r++) {
for (unsigned c = 0; c < nc; c++) {
mv[r*nc + c] = 0;
}
}
}
void
echelon_form_bareiss_skip_statfactors(matrix &m)
{
unsigned nr = m.rows();
unsigned nc = m.cols();
exvector &mv = ((matrix_hack*)&m)->get_m();
exvector extraf(nr);
for (unsigned r = 0; r < nr; r++) { extraf[r] = 1; }
ex z = 1;
unsigned r0 = 0;
for (unsigned c0 = 0; (c0 < nc) && (r0 < nr - 1); c0++) {
unsigned pivot = r0;
// No normalization before is_zero() here, because
// we maintain the matrix normalized throughout the
// algorithm.
while ((pivot < nr) && mv[pivot*nc + c0].is_zero()) { pivot++; }
if (pivot == nr) {
// The whole column below r0:c0 is zero, we can skip
// it.
continue;
}
if (pivot > r0) {
// Found a non-zero row somewhere below r0; let's
// swap it in.
for (unsigned c = c0; c < nc; c++) {
mv[pivot*nc + c].swap(mv[r0*nc + c]);
}
extraf[pivot].swap(extraf[r0]);
}
ex a = mv[r0*nc + c0];
for (unsigned r = r0 + 1; r < nr; r++) {
ex b = mv[r*nc + c0];
if (b.is_zero()) {
// This branch is the same as the next one, but
// is is simplified due to b being 0.
ex f;
assert(divide(a*extraf[r0]*extraf[r], z, f));
extraf[r] = f;
} else {
ex g = gcd(a*extraf[r0], b*extraf[r]);
ex f;
gcd(extraf[r0]*g, z, &f);
f = normal(f);
ex d;
assert(divide(z*f, extraf[r0]*extraf[r], d));
for (unsigned c = c0 + 1; c < nc; c++) {
ex v;
assert(divide(a*mv[r*nc + c] - b*mv[r0*nc + c], d, v));
// Since v is a polynomial, the call to normal()
// here does not involve GCD, but rather just
// expands sums while preserving common factors,
// which is nomally preferable to full expand().
m[r*nc + c] = normal(v);
}
mv[r*nc + c0] = 0;
extraf[r] = f;
}
}
z = a*extraf[r0];
r0++;
}
// Zero out the remaining rows (just in case).
for (unsigned r = r0 + 1; r < nr; r++) {
for (unsigned c = 0; c < nc; c++) {
mv[r*nc + c] = 0;
}
}
// At this point mv[r*nc + c]*extraf[r] is precisely what
// plain Bareiss algorithm would have computed; we can now
// discard this extra factor freely.
}
void
save_matrix(ostream &f, const matrix &m)
{
f << "{";
for (unsigned i = 0; i < m.rows(); i++) {
if (i != 0) f << ",\n";
f << "{";
for (unsigned j = 0; j < m.cols(); j++) {
if (j != 0) f << ", ";
f << m(i, j);
}
f << "}";
}
f << "}";
}
int
main()
{
symbol a11("a11"), a12("a12"), a13("a13"), a14("a14"), a15("a15");
symbol a21("a21"), a22("a22"), a23("a23"), a24("a24"), a25("a25");
symbol a31("a31"), a32("a32"), a33("a33"), a34("a34"), a35("a35");
symbol a41("a41"), a42("a42"), a43("a43"), a44("a44"), a45("a45");
symbol a51("a51"), a52("a52"), a53("a53"), a54("a54"), a55("a55");
matrix m = {
{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}
};
cout << "original matrix:" << endl;
save_matrix(cout, ex_to<matrix>(factor(m)));
cout << endl << endl;
matrix m1 = m;
echelon_form_bareiss_original(m1);
cout << "bareiss original:" << endl;
save_matrix(cout, ex_to<matrix>(factor(m1)));
cout << endl << endl;
matrix m2 = m;
echelon_form_bareiss_skip_zrows(m2);
cout << "bareiss/skip zrows:" << endl;
save_matrix(cout, ex_to<matrix>(factor(m2)));
cout << endl << endl;
matrix m3 = m;
echelon_form_bareiss_skip_statfactors(m3);
cout << "bareiss/skip statfactors:" << endl;
save_matrix(cout, ex_to<matrix>(factor(m3)));
cout << endl;
}
More information about the GiNaC-list
mailing list