Skip to content

Commit 5351d0c

Browse files
The determinant is computed by fraction-free elimination where it can be (#1129)
`Matrix.Determinant` expanded by Laplace, which is O(n!). For a fully symbolic matrix that is optimal -- the determinant genuinely has n! terms and no algorithm returns it smaller in expanded form. For a numeric one it is pure waste: the answer is a single number and O(n^3) work suffices. Bareiss' fraction-free elimination now runs wherever the entries are polynomials over the rationals, and Laplace answers everything else. What decides it is not the size but whether the entries can be read, settled per matrix by trying -- a 12x12 with two variables and a 4x4 with sixteen are different problems and n alone does not tell them apart. Numeric, both arms built from source on one machine: 8x8 from 382 ms to under 1 ms, 10x10 from 14 415 ms to 2 ms, and 11x11, 12x12, 20x20 and 30x30 from not returning at all to 2, 3, 11 and 22 ms. **The elimination was already here.** `PolynomialResultant.SylvesterDeterminant` is a Bareiss elimination over `MultivariatePolynomial`, written for the Sylvester matrix and carrying everything this needs: exact division, row swaps, sign tracking, a work budget and a cancellation check. It is split into `FractionFreeDeterminant` and shared rather than written twice, which is worth saying out loud: a sign convention or an off-by-one in an elimination is a wrong answer that looks entirely plausible, and one implementation exercised by two callers is tested by both. **No condition is introduced.** The issue is right that this is the thing to check rather than argue: exactness is a fact about the ring, and `Entity` division would not notice it -- `(a * b) / a` stays written as a quotient. The arithmetic here never touches `Entity` division. It happens in `MultivariatePolynomial`, which has no quotients to leave behind, and `DivideExact` returns null rather than a remainder, so a division that does not come out stops the elimination and sends the caller to Laplace. The answer is a polynomial in the entries, or it is Laplace's. Declined, and answered by Laplace exactly as before: an entry that is not a polynomial over the rationals, a matrix in more than eight indeterminates, and a matrix mentioning `e` or `pi`, since a constant is a value rather than an indeterminate and this ring cannot hold one. The printed form of a symbolic determinant changes, because an elimination produces an expanded polynomial where Laplace produces a nested one. Recorded in BREAKING-CHANGES.md with the shapes measured on both arms; the expression is no larger in any case measured, and `x * (x ^ 2 + -1) + -x` becomes `x ^ 3 - 2 * x`. The two algorithms are compared on 300 generated matrices where both apply, as a difference simplified to zero rather than as trees -- they group the same polynomial differently and that is not a disagreement. No disagreements. The comparison also asserts that it read more than a third of what it generated, so a change making the elimination decline everything fails rather than passes vacuously. Part of #999.
1 parent 7a68049 commit 5351d0c

7 files changed

Lines changed: 539 additions & 59 deletions

File tree

‎BREAKING-CHANGES.md‎

Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,8 @@ read first.
2828
| **Silent** | `MathS.Abs("x").WithCodomain(Domain.Any).Stringize()`, and every node widened to `Any` from a narrower default | `abs(x)` — reads back as `Real`, losing the widening | `domain(abs(x), Any)` |
2929
| **Silent** | `"domain(1/2, CC)".ToEntity()`, and every quotient of two integer literals annotated with the codomain its node type does not default to | `1/2` — the annotation dropped, equal to the unannotated literal | `1/2` carrying `Codomain = Complex`, which prints and reads back as `domain(1/2, CC)` |
3030
| | `"a in (a / 3; 3a)".ToEntity().Simplify()`, and every denominator but 2 | `a in (a / 3; 3 * a)` — left as written | `a > 0` |
31+
| | `MathS.Matrix(...).Determinant` on a symbolic matrix of polynomials | `a * d + -b * c`, Laplace's nested expansion | `a * d - b * c`, expanded — the same value, no larger |
32+
| | the same on a numeric matrix past 10x10 | did not return | 11x11 in 2 ms, 30x30 in 22 ms |
3133
| | `"a in (a / 2; 0)".ToEntity().Simplify()`, and every interval demanding both signs | left as written | `False provided a in RR` |
3234
| | `new Entity[0].SumAll()`, and `Sumf.Sum` on an empty list | `AngouriBugException: At least 1 child required` | `0` |
3335
| | `new Entity[0].MultiplyAll()`, and `Mulf.Multiply` on an empty list | `AngouriBugException` | `1` |
@@ -848,6 +850,57 @@ together because the first is what makes the second free: the guard alone cost t
848850
that were being answered by coincidence, and the collection restores them along with the rest of the
849851
family ([#1056](https://github.com/asc-community/AngouriMath/issues/1056)).
850852

853+
### The determinant is computed by fraction-free elimination where it can be
854+
855+
`Matrix.Determinant` expanded by Laplace, which is `O(n!)`. For a fully symbolic matrix that is
856+
optimal — the determinant genuinely has `n!` terms, and no algorithm returns it smaller in expanded
857+
form. For a numeric one it is pure waste: the answer is a single number and `O(n^3)` work suffices.
858+
859+
Bareiss' fraction-free elimination now runs wherever the entries are polynomials over the rationals,
860+
and Laplace answers everything else. What decides it is not the size but whether the entries can be
861+
read, settled per matrix by trying.
862+
863+
**The ceiling this removes**, both arms built from source on one machine:
864+
865+
| | before | now |
866+
|---|---|---|
867+
| numeric 8×8 | 382 ms | under 1 ms |
868+
| numeric 10×10 | 14 415 ms | 2 ms |
869+
| numeric 11×11 | did not return in four minutes | 2 ms |
870+
| numeric 12×12 | did not return | 3 ms |
871+
| numeric 20×20 | did not return | 11 ms |
872+
| numeric 30×30 | did not return | 22 ms |
873+
874+
**The printed form of a symbolic determinant changes**, because an elimination produces an expanded
875+
polynomial where Laplace produces a nested expansion. The value is the same and the expression is no
876+
larger in any case measured:
877+
878+
| | before | now |
879+
|---|---|---|
880+
| `"[[a, b], [c, d]]"` | `a * d + -b * c` | `a * d - b * c` |
881+
| `"[[x, 1, 0], [1, x, 1], [0, 1, x]]"` | `x * (x ^ 2 + -1) + -x` | `x ^ 3 - 2 * x` |
882+
| `"[[a, b, 1], [c, d, 2], [1, 2, 3]]"` | `a * (d * 3 + -4) + -b * (c * 3 + -2) + c * 2 + -d` | `3 * a * d - 4 * a - 3 * b * c + 2 * b + 2 * c - d` |
883+
| `"[[x, 1], [1, x]]"` | `x ^ 2 + -1` | unchanged |
884+
| `"[[1/2, 1/3], [1/4, 1/5]]"` | `1/60` | unchanged |
885+
886+
**No condition is introduced, and that is the point.** An ordinary Gaussian elimination leaves its
887+
pivots as literal divisions, so its answer is undefined wherever a pivot vanishes — at points where
888+
the determinant is perfectly well defined. That was
889+
[#992](https://github.com/asc-community/AngouriMath/issues/992), and it is why Laplace was chosen.
890+
Bareiss divides as well, but each division is by the *previous* pivot and is exact: the quotient is a
891+
determinant of a minor, so it is back in the ring. Here it is exact **and checked** — the arithmetic
892+
happens in `MultivariatePolynomial`, which has no quotients to leave behind, and a division that does
893+
not come out returns null and sends the caller to Laplace.
894+
895+
**What is declined**, and answered by Laplace exactly as before: an entry that is not a polynomial
896+
over the rationals (`sin(x)`, `1 / x`, `2 ^ x`), a matrix in more than eight indeterminates, and a
897+
matrix mentioning `e` or `pi` — a constant is a value rather than an indeterminate, and this ring
898+
cannot hold one.
899+
900+
The two algorithms were compared on 300 generated matrices where both apply, as a difference
901+
simplified to zero rather than as trees, with **no disagreements**
902+
([#999](https://github.com/asc-community/AngouriMath/issues/999)).
903+
851904
### A fold over an empty sequence answers instead of throwing
852905

853906
A fold over a monoid has an identity, so the empty sum is `0` and the empty product is `1` — which

‎Sources/.editorconfig‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -51,6 +51,15 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed
5151
[AngouriMath/Functions/Algebra/Polynomials/KroneckerFactorization.cs]
5252
file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n
5353

54+
[AngouriMath/Functions/Algebra/Polynomials/FractionFreeDeterminant.cs]
55+
file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n
56+
57+
[AngouriMath/Functions/Algebra/Polynomials/PolynomialDeterminant.cs]
58+
file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n
59+
60+
[Tests/UnitTests/Algebra/Polynomials/FractionFreeDeterminantTest.cs]
61+
file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n
62+
5463
[AngouriMath/Functions/Algebra/Groebner/*.cs]
5564
file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n
5665

‎Sources/AngouriMath/Core/Entity/Omni/Entity.Matrix.cs‎

Lines changed: 25 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -270,23 +270,37 @@ public Entity AsScalar()
270270
this);
271271
private LazyPropertyA<Matrix> t;
272272

273-
// We do not need to use Gaussian elimination here
274-
// since we anyway get N! memory use.
275-
// Gaussian elimination is also not merely no cheaper, it is wrong here:
276-
// it leaves the pivots as literal divisions, so the expression it returns
277-
// is undefined wherever a pivot vanishes -- at points where the
278-
// determinant itself is perfectly well defined. The determinant of a
279-
// matrix over a commutative ring is a polynomial in its entries, and
280-
// Laplace expansion never divides, so there is nothing to exclude.
273+
// Two algorithms, and which one runs is decided by the entries rather than by the
274+
// size. Bareiss' fraction-free elimination is O(n^3) and Laplace expansion is
275+
// O(n!), so Bareiss is tried first -- but it needs entries it can read as
276+
// polynomials over the rationals, and it says so by declining rather than by
277+
// guessing. https://github.com/asc-community/AngouriMath/issues/999
278+
//
279+
// *Ordinary* Gaussian elimination is not merely no cheaper than Laplace, it is
280+
// wrong here: it leaves the pivots as literal divisions, so the expression it
281+
// returns is undefined wherever a pivot vanishes -- at points where the
282+
// determinant itself is perfectly well defined. The determinant of a matrix over a
283+
// commutative ring is a polynomial in its entries, and neither Laplace nor Bareiss
284+
// ever leaves a quotient, so there is nothing to exclude either way.
281285
// https://github.com/asc-community/AngouriMath/issues/992
282286
/// <summary>
283-
/// Finds the symbolical determinant via Laplace's method
287+
/// The determinant, by fraction-free elimination where the entries are polynomials
288+
/// over the rationals and by Laplace expansion otherwise.
284289
/// </summary>
290+
/// <remarks>
291+
/// Both are exact and neither divides in a way that could exclude a point, so the
292+
/// two agree wherever both apply; which one ran is a fact about cost, not about the
293+
/// answer. Laplace is optimal for a fully symbolic matrix, whose determinant
294+
/// genuinely has n! terms.
295+
/// </remarks>
285296
public Entity? Determinant => determinant.GetValue(
286-
static @this =>
297+
static @this =>
287298
{
288-
if (!@this.IsSquare)
299+
if (!@this.IsSquare)
289300
return null;
301+
if (Functions.PolynomialDeterminant.Of(@this.RowCount, (r, c) => @this[r, c])
302+
is { } byElimination)
303+
return byElimination.InnerSimplified;
290304
return @this.InnerMatrix.DeterminantLaplace().InnerSimplified;
291305
},
292306
this
Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,127 @@
1+
//
2+
// Copyright (c) 2019-2026 Angouri.
3+
// AngouriMath is licensed under MIT.
4+
// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.
5+
// Website: https://am.angouri.org.
6+
//
7+
8+
using AngouriMath.Core.Multithreading;
9+
using PeterO.Numbers;
10+
11+
namespace AngouriMath.Functions
12+
{
13+
/// <summary>
14+
/// Bareiss' fraction-free elimination: the determinant of a square matrix over an
15+
/// integral domain, in <c>O(n^3)</c> ring operations and without ever leaving the ring.
16+
/// </summary>
17+
/// <remarks>
18+
/// <para>
19+
/// Ordinary Gaussian elimination divides by the pivot, so over a ring of polynomials it
20+
/// produces quotients — and an expression built from them is undefined wherever a pivot
21+
/// vanishes, at points where the determinant itself is perfectly well defined. That is
22+
/// what <a href="https://github.com/asc-community/AngouriMath/issues/992">#992</a> was
23+
/// about, and why <see cref="Entity.Matrix.Determinant"/> used Laplace expansion rather
24+
/// than elimination.
25+
/// </para>
26+
/// <para>
27+
/// Bareiss divides too, and the divisions are the point: each entry is divided by the
28+
/// <em>previous</em> pivot, and that division is <b>exact</b> — the quotient is a
29+
/// determinant of a minor and therefore back in the ring. So the intermediate entries do
30+
/// not swell the way Gaussian elimination's do, and nothing is ever left as a quotient to
31+
/// exclude a point.
32+
/// </para>
33+
/// <para>
34+
/// <b>Exactness is a fact about the ring, and is checked here anyway.</b>
35+
/// <see cref="MultivariatePolynomial.DivideExact(MultivariatePolynomial, int)"/> returns
36+
/// null rather than a remainder, so a division that does not come out — because the term
37+
/// ceiling was reached, or because an assumption above is wrong — stops the elimination
38+
/// instead of producing an answer that is quietly not the determinant. Every caller reads
39+
/// null as "use another method", never as "there is no determinant".
40+
/// </para>
41+
/// <para>
42+
/// Split out of <see cref="PolynomialResultant"/>, where it was written for the Sylvester
43+
/// matrix. Shared rather than written twice, deliberately: a sign convention or an
44+
/// off-by-one in an elimination is a wrong answer that looks entirely plausible, and one
45+
/// implementation exercised by two callers is tested by both.
46+
/// <a href="https://github.com/asc-community/AngouriMath/issues/999">#999</a>
47+
/// </para>
48+
/// </remarks>
49+
internal static class FractionFreeDeterminant
50+
{
51+
[ConstantField] private static readonly ERational MinusOne = ERational.One.Negate();
52+
53+
/// <summary>
54+
/// The determinant of <paramref name="matrix"/>, or <see langword="null"/> where the
55+
/// elimination could not finish within <paramref name="maxWork"/> or a division came
56+
/// back inexact.
57+
/// </summary>
58+
/// <param name="matrix">
59+
/// Square, and <b>modified in place</b> — the caller owns it and must not need it
60+
/// afterwards. Elimination is destructive, and copying here would hide the cost of the
61+
/// one thing this exists to make cheap.
62+
/// </param>
63+
/// <param name="variableCount">
64+
/// How many variables the entries are over, needed to build the zero and the one of the
65+
/// same ring. Not derivable from the entries, since a matrix may be all constants.
66+
/// </param>
67+
/// <param name="maxWork">
68+
/// A ceiling on term-pair multiplications, charged <em>before</em> each multiplication
69+
/// rather than after, so that the budget cannot be overshot by the one step that was
70+
/// going to be the most expensive of them.
71+
/// </param>
72+
internal static MultivariatePolynomial? Of(
73+
MultivariatePolynomial[][] matrix, int variableCount, long maxWork)
74+
{
75+
var size = matrix.Length;
76+
if (size == 0)
77+
return MultivariatePolynomial.One(variableCount);
78+
79+
var zero = MultivariatePolynomial.Zero(variableCount);
80+
var negated = false;
81+
var work = 0L;
82+
var previous = MultivariatePolynomial.One(variableCount);
83+
84+
for (var pivot = 0; pivot + 1 < size; pivot++)
85+
{
86+
MultithreadingFunctional.ExitIfCancelled();
87+
if (matrix[pivot][pivot].IsZero)
88+
{
89+
var replacement = -1;
90+
for (var row = pivot + 1; row < size && replacement < 0; row++)
91+
if (!matrix[row][pivot].IsZero)
92+
replacement = row;
93+
// Nothing below the pivot to bring up means the column is a combination
94+
// of the ones before it, so the matrix is singular and its determinant
95+
// is zero. That is an answer, not a failure.
96+
if (replacement < 0)
97+
return zero;
98+
(matrix[pivot], matrix[replacement]) = (matrix[replacement], matrix[pivot]);
99+
negated = !negated;
100+
}
101+
var head = matrix[pivot][pivot];
102+
for (var row = pivot + 1; row < size; row++)
103+
{
104+
MultithreadingFunctional.ExitIfCancelled();
105+
var leading = matrix[row][pivot];
106+
for (var column = pivot + 1; column < size; column++)
107+
{
108+
work += (long)head.TermCount * matrix[row][column].TermCount
109+
+ (long)leading.TermCount * matrix[pivot][column].TermCount;
110+
if (work > maxWork)
111+
return null;
112+
if (head.Multiply(matrix[row][column]) is not { } kept
113+
|| leading.Multiply(matrix[pivot][column]) is not { } removed
114+
|| kept.Subtract(removed).DivideExact(previous) is not { } reduced)
115+
return null;
116+
matrix[row][column] = reduced;
117+
}
118+
matrix[row][pivot] = zero;
119+
}
120+
previous = head;
121+
}
122+
123+
var determinant = matrix[size - 1][size - 1];
124+
return negated ? determinant.ScaleBy(MinusOne) : determinant;
125+
}
126+
}
127+
}

0 commit comments

Comments
 (0)