half of rowEchelon

[Ralf Hemmecke]

Hello,

the result of rowEchelon is nice, but a bit too much for what I need.
Is there a function that just kills the elements below the diagonal and
doesn't care about the upper triangular part.

So I would be happy with rowEchelon with a replacement of the loop
coming at the end

https://github.com/fricas/fricas/blob/master/src/algebra/matfuns.spad#L119

-- perform row operations so that jth column has only one 1
for k in minR..maxR repeat
if k ~= i and not((pp := qelt(x, k, j)) = 0) then
row_op(x, i, j + 1, k, maxC, pp)
qsetelt!(x, k, j, 0)

by

for k in i+1..maxR repeat
if not((pp := qelt(x, k, j)) = 0) then
row_op(x, i, j + 1, k, maxC, pp)
qsetelt!(x, k, j, 0)

Of course, such a function needs a new name, and the above covers only
the Field case not EuclideanDomain.

Any ideas about what to do?

Ralf



(298) -> mat := matrix [[1,1,1,1],[1,5,3,4],[1,4,9,16],[2,3,5,7]]

+1 1 1 1 +
| |
|1 5 3 4 |
(298) | |
|1 4 9 16|
| |
+2 3 5 7 +
(300) -> rowEchelon mat

+1 0 0 13+
| |
|0 1 0 5 |
(300) | |
|0 0 1 0 |
| |
+0 0 0 17+
Type: Matrix(Integer)

(301) -> rowEchelon(mat::Matrix(Fraction Integer))

+1 0 0 0+
| |
|0 1 0 0|
(301) | |
|0 0 1 0|
| |
+0 0 0 1+
Type: Matrix(Fraction(Integer))

[Waldek Hebisch]

On Fri, May 12, 2023 at 02:03:15PM +0200, Ralf Hemmecke wrote:
> Hello,
>
> the result of rowEchelon is nice, but a bit too much for what I need.
> Is there a function that just kills the elements below the diagonal and
> doesn't care about the upper triangular part.
>
> So I would be happy with rowEchelon with a replacement of the loop coming at
> the end
>
> https://github.com/fricas/fricas/blob/master/src/algebra/matfuns.spad#L119
>
> -- perform row operations so that jth column has only one 1
> for k in minR..maxR repeat
> if k ~= i and not((pp := qelt(x, k, j)) = 0) then
> row_op(x, i, j + 1, k, maxC, pp)
> qsetelt!(x, k, j, 0)
>
> by
>
> for k in i+1..maxR repeat
> if not((pp := qelt(x, k, j)) = 0) then
> row_op(x, i, j + 1, k, maxC, pp)
> qsetelt!(x, k, j, 0)
>
> Of course, such a function needs a new name, and the above covers only the
> Field case not EuclideanDomain.
>
> Any ideas about what to do?

I am not sure if understand what you actually want. In general
"row echelon" form is quite non-unique. To make it better determined
we require that all elements above diagonal are reduced as
much as possible, which in case of working over a field means
reduced to zero in columns containing pivots

So do you want "row echelon" form but are worried about cost of
zeroing elements above diagonal? Or do you want to preserve
information in upper part?

You probably should look at LU decomposition ('LUDecomp' in
package LUDecomposition). 'LUDecomp' is oriented towards
equation solving and produces more results than you apparently
want, but in a sense performs minimal number of operations
needed to turn elements below diagonal into zeros.
'LUDecomp' has restrictions (I think that it wants nonzero
pivot in each column), but there is generalization without
such restriction (I have uncommited version which I developed
for unfinished modular solver for linear systems).

--
Waldek Hebisch

[Ralf Hemmecke]

> So do you want "row echelon" form but are worried about cost of
> zeroing elements above diagonal? Or do you want to preserve
> information in upper part?

Well, my matrix actually comes from polynomial coefficients of several
polynomials. I want them to be of different degree eventually and want
to record how these different-degree polynomials came from the original
ones. Roughly speaking I want the "representation" of the new in terms
of the old as "small" as possible. Any tail-reduction would add a new
term to the representation, so I would like to avoid this.

> You probably should look at LU decomposition ('LUDecomp' in
> package LUDecomposition).

I've meanwhile found fractionFreeGauss! in MatrixLinearAlgebraFunctions.
That seems to do what I want. It even looks as if it divides out factors
that were multiplied in a previous step (the "exquo b" part) here.

https://github.com/fricas/fricas/blob/master/src/algebra/matfuns.spad#L623
https://github.com/fricas/fricas/blob/master/src/algebra/matfuns.spad#L628

The only thing I do not yet understand is what the ans thing does?
Why would I have to negate the last row in case of odd row swaps?

Ralf

[Waldek Hebisch]

On Fri, May 12, 2023 at 10:36:48PM +0200, Ralf Hemmecke wrote:
> > So do you want "row echelon" form but are worried about cost of
> > zeroing elements above diagonal? Or do you want to preserve
> > information in upper part?
>
> Well, my matrix actually comes from polynomial coefficients of several
> polynomials. I want them to be of different degree eventually and want to
> record how these different-degree polynomials came from the original ones.
> Roughly speaking I want the "representation" of the new in terms of the old
> as "small" as possible. Any tail-reduction would add a new term to the
> representation, so I would like to avoid this.

Let me remark that what you write above looks a bit similar
to Hermite-Pade problem and we have specialized solvers for
Hermite-Pade problem. If you have different problem, then
Hermite-Pade probably will not help.

> > You probably should look at LU decomposition ('LUDecomp' in
> > package LUDecomposition).
>
> I've meanwhile found fractionFreeGauss! in MatrixLinearAlgebraFunctions.
> That seems to do what I want.

OK.

> It even looks as if it divides out factors
> that were multiplied in a previous step (the "exquo b" part) here.
>
> https://github.com/fricas/fricas/blob/master/src/algebra/matfuns.spad#L623
> https://github.com/fricas/fricas/blob/master/src/algebra/matfuns.spad#L628

That looks like what is usually called Bareiss elimination.
Without such step it would be unusable for larger problems
due to exponential complexity (exponential growth of size of
coefficients).

> The only thing I do not yet understand is what the ans thing does?
> Why would I have to negate the last row in case of odd row swaps?

'ans' is to preserve sign of determinanat. Row swap multiplies sign
of determinant by -1, so we need to compensate this.

--
Waldek Hebisch

[Ralf Hemmecke]

>> I've meanwhile found fractionFreeGauss! in MatrixLinearAlgebraFunctions.
>> That seems to do what I want.

> That looks like what is usually called Bareiss elimination.

Yes, exactly. Unfortunately, users might know about Bareiss, but to find
fractionFreeGaus! is hard, since nowhere in the FriCAS code "bareiss" is
mentioned. I think we should add it to the documentation of
fractionFreeGaus! that this is the Bareiss algorithm.

> 'ans' is to preserve sign of determinanat. Row swap multiplies sign
> of determinant by -1, so we need to compensate this.

Yes, now it's clear.

Would you mind if I propose a patch with some more hints in the code?
(see pull request at github)

https://github.com/fricas/fricas/pull/123


Ralf

[Waldek Hebisch]

OK

--
Waldek Hebisch