# How to exactly solve large overdetermined systems of linear equations
**URL:** https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710
**Category:** Numerics
**Created:** [September 3, 2022, 12:21am UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710 "2022-09-03T00:21:36Z")
**Posts on this page:** 17
**Page:** 1
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [September 3, 2022, 12:21am UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/1 "2022-09-03T00:21:36Z")
Suppose I have, e.g., 250 linear equations with `Rational{BigInt}` coefficients, with 150 variables. If the system is inconsistent I want to know about it, otherwise I want the _exact_ solution.
Is there an easier way than writing my own Gaussian elimination?
This won’t run often, so I’m not concerned with long runtimes.
---
### Author:  [@jd-foster](https://discourse.julialang.org/u/jd-foster)
#### Post date: [September 3, 2022, 12:30am UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/2 "2022-09-03T00:30:22Z")
Have a look at Nemo.jl; it does exact linear algebra.
---
### Author:  [@thofma](https://discourse.julialang.org/u/thofma)
#### Post date: [September 5, 2022, 12:04pm UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/3 "2022-09-05T12:04:30Z")
Just for sake of completeness, either using AbstractAlgebra or (more efficiently) Nemo, it can be done as follows:
```julia
julia> using AbstractAlgebra # will also work with using Nemo;
julia> A = QQ[1 0; 3 0; 5 0]; # create an AbstractAlgebra matrix
julia> v = QQ[1; 2; 3]; # create a 3x1 matrix
julia> can_solve_with_solution(A, v) # solve Ax = v
(false, [1; 0])
julia> A = QQ[1 2; 3 4; 5 6];
julia> can_solve_with_solution(A, v) # solve Ax = v
(true, [0; 1//2])
julia> fl, x = can_solve_with_solution(A, v) # solve Ax = v
(true, [0; 1//2])
julia> A*x == v
true
```
One can also do `can_solve_with_solution(A, v; side = :left)` to solve `xA = v`.
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [September 5, 2022, 12:28pm UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/4 "2022-09-05T12:28:57Z")
TBH I didn’t even realize that `can_solve_with_solution` existed, instead i used `AbstractAlgebra.rref` 😅
---
### Author:  [@gustave\_Robichon](https://discourse.julialang.org/u/gustave_Robichon)
#### Post date: [March 5, 2024, 4:43pm UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/5 "2024-03-05T16:43:20Z")
hello thofma,
this interests me very much. does a QQMatrix has some elements of the type BigInt ? I want to solve a very big system A.X=Y with Rational{BigInt}. Thanks !
And, is there a fast way to convert a Matrix{Rational{BigInt}} to a QQMatrix, and back ? I could copy it term by term, eventually …
thanks a lot 🙂
---
### Author:  [@thofma](https://discourse.julialang.org/u/thofma)
#### Post date: [March 5, 2024, 9:53pm UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/6 "2024-03-05T21:53:12Z")
Yes, there is a fast way. Here is an example, where one solves `A*X = Y` where `X` and `Y` are of type `Vector`.
```julia-repl
julia> AA = Rational{BigInt}[1 2; 3 4]
2×2 Matrix{Rational{BigInt}}:
1 2
3 4
julia> A = matrix(QQ, AA)
[1 2]
[3 4]
julia> YY = Rational{BigInt}[5, 6]
2-element Vector{Rational{BigInt}}:
5
6
julia> Y = QQ.(YY)
2-element Vector{QQFieldElem}:
5
6
julia> fl, X = can_solve_with_solution(A, Y; side = :right) # solves A * X = Y for X
(true, QQFieldElem[-4, 9//2])
julia> Rational{BigInt}.(X)
2-element Vector{Rational{BigInt}}:
-4
9//2
```
Here is an example, where `X` and `Y` are themselves of type `Matrix`:
```julia-repl
julia> YY = Rational{BigInt}[5 6; 7 8]
2×2 Matrix{Rational{BigInt}}:
5 6
7 8
julia> Y = matrix(QQ, YY)
[5 6]
[7 8]
julia> fl, X = can_solve_with_solution(A, Y; side = :right) # solves A * X = Y for X
(true, [-3 -4; 4 5])
julia> XX = Matrix{Rational{BigInt}}(X)
2×2 Matrix{Rational{BigInt}}:
-3 -4
4 5
```
Depending on how what you mean with “very big”, it might take a while to solve your system, but give it a try!
---
### Author:  [@gustave\_Robichon](https://discourse.julialang.org/u/gustave_Robichon)
#### Post date: [March 7, 2024, 9:29am UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/7 "2024-03-07T09:29:40Z")
Ok thanks, everything works well 🙂 I am new to Julia, and I don’t understand “constructors” well !
“Big Matrix” means only that I might need some BigInt, in any case the coefficients explode 🙂
Thanks a lot !
Gustave
---
### Author:  [@gustave\_Robichon](https://discourse.julialang.org/u/gustave_Robichon)
#### Post date: [March 7, 2024, 1:20pm UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/8 "2024-03-07T13:20:22Z")
Hello 🙂
I actually I should use Matrix{Complex{Rational{Int}}} . Is there way to do so with Nemo ? Or maybe with AbstractAlgebra ? Thanks ! 🙂
finally I did the code with Complex{Float64}. But I like Algebra, so that your answer interests me !!!
Bye 🙂
---
### Author:  [@thofma](https://discourse.julialang.org/u/thofma)
#### Post date: [March 16, 2024, 10:15pm UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/9 "2024-03-16T22:15:27Z")
Sorry for the late reply. Yes, this is also possible:
```julia-repl
julia> QQi = Nemo.QQiField()
Gaussian rational field
julia> AA = Complex{Rational{BigInt}}[1 2im; 3 4im]
2×2 Matrix{Complex{Rational{BigInt}}}:
1//1+0//1*im 0//1+2//1*im
3//1+0//1*im 0//1+4//1*im
julia> A = matrix(QQi, AA)
[1 2*im]
[3 4*im]
julia> YY = Complex{Rational{BigInt}}[5im, 6]
2-element Vector{Complex{Rational{BigInt}}}:
0//1 + 5//1*im
6//1 + 0//1*im
julia> Y = QQi.(YY)
2-element Vector{Nemo.QQiFieldElem}:
5*im
6
julia> fl, X = can_solve_with_solution(A, Y; side = :right) # solves A * X = Y for X
(true, Nemo.QQiFieldElem[6 - 10*im, 15//2 + 3*im])
julia> XX = [Complex{Rational{BigInt}}(real(y), imag(y)) for y in X]
2-element Vector{Complex{Rational{BigInt}}}:
6//1 - 10//1*im
15//2 + 3//1*im
```
---
### Author:  [@gustave\_Robichon](https://discourse.julialang.org/u/gustave_Robichon)
#### Post date: [March 17, 2024, 7:45am UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/10 "2024-03-17T07:45:08Z")
OK thanks, it’s perfect ! I gonna try it tomorrow 🙂
Have a nice day 🙂
Gustave
---
### Author:  [@gustave\_Robichon](https://discourse.julialang.org/u/gustave_Robichon)
#### Post date: [March 18, 2024, 9:41am UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/11 "2024-03-18T09:41:31Z")
It works perfectly well, thank you !!!
---
### Author:  [@gustave\_Robichon](https://discourse.julialang.org/u/gustave_Robichon)
#### Post date: [January 2, 2025, 3:31pm UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/12 "2025-01-02T15:31:25Z")
hello, once again 🙂
I want to solve AX=Y, for A a matrix and Y a vector, for some Rational{Polynomial{Complex{Rational{BigInt}}}}, is it possible ? I found everything in the documentation of AbstractAlgebra, but the Complex{ } part. Any help ?
I’ve also tried with A \ Y, where the awkward type is :
> RationalPoly{Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Complex{Rational{BigInt}}},
> Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Complex{Rational{BigInt}}}}}
generated by “@polyvar eta\_var”
but strangely it works for a small matrix but not for a big one. I think it changes the algorithm when the size increases, or I didn’t understand what happens there.
It’s frustrating because this is a simple algorithm, but I can’t find it anywhere ☹ Maybe I could write it by myself.
Thanks for any help !
Gustave
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [January 2, 2025, 4:28pm UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/13 "2025-01-02T16:28:52Z")
> [@gustave\_Robichon](#):
>
> I want to solve AX=Y, for A a matrix and Y a vector, for some Rational{Polynomial{Complex{Rational{BigInt}}}}, is it possible ? I found everything in the documentation of AbstractAlgebra, but the Complex{ } part. Any help ?
```julia-repl
julia> using Nemo
Welcome to Nemo version 0.48.0
Nemo comes with absolutely no warranty whatsoever
julia> QQi = Nemo.GaussianRationals()
Gaussian rational field
julia> QQi(3, 4)
3 + 4*im
julia> QQi(3 + 4im)
3 + 4*im
```
> [@gustave\_Robichon](#):
>
> but strangely it works for a small matrix but not for a big one. I think it changes the algorithm when the size increases, or I didn’t understand what happens there.
Sounds like a bug. Can you provide a reproducer, so the bug could be reported?
---
### Author:  [@gustave\_Robichon](https://discourse.julialang.org/u/gustave_Robichon)
#### Post date: [January 3, 2025, 10:44am UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/14 "2025-01-03T10:44:42Z")
Thank you !
For the bug, I describe what I use, what works and what does not.
```julia
using DynamicPolynomials
@polyvar eta_var
eta = (BigInt(0)//1 + (1//1 + 1im*0//1) *eta_var) / (BigInt(1)//1 +(0//1+0//1 * 1im) *eta_var)
```
the goal of using eta is to have a type of a rational{Polynomial{Rational{Complex}}
the type of eta is :
```julia
RationalPoly{Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Complex{Rational{BigInt}}}, Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Complex{Rational{BigInt}}}}
```
```julia
A = [eta 2*eta ; 0 eta * eta]
YA = [(1+0 * eta) (1+0 * eta)] '
```
A \ YA works
```julia
B # Is a matrix
21×10 Matrix{RationalPoly{Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Complex{Rational{BigInt}}}, Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Complex{Rational{BigInt}}}}}:
YB # Is a vector
21-element Vector{RationalPoly{Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Complex{Rational{BigInt}}}, Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Complex{Rational{BigInt}}}}}:
```
```julia
B[1,2] = ((-1//8 + 0//1*im) eta_var^2) / ((1//1 + 0//1*im))
B[3,2] = ((0//1 + 1//4*im)*eta_var) / ((1//1 + 0//1*im))
B[4,3] = ((-1//4 + 0//1*im)*eta_var^2) / ((1//1 + 0//1*im))
B[5,3] = ((0//1 - 1//4*im)*eta_var) / ((1//1 + 0//1*im))
B[6,3] = ((0//1 + 1//4*im)*eta_var) / ((1//1 + 0//1*im))
B[7,4] = ((-1//8 + 0//1*im)*eta_var^2) / ((1//1 + 0//1*im))
B[8,4] = ((0//1 + 1//4*im)*eta_var) / ((1//1 + 0//1*im))
B[10,1] = ((1//1 + 0//1*im)) / ((1//1 + 0//1*im))
B[11,6] = ((-1//4 + 0//1*im)) / ((1//1 + 0//1*im))
B[12,7] = ((-1//4 + 0//1*im)) / ((1//1 + 0//1*im))
B[13,8] = ((-1//8 + 0//1*im)*eta_var^2) / ((1//1 + 0//1*im))
B[14,8] = ((0//1 - 1//4*im)*eta_var) / ((1//1 + 0//1*im))
B[16,9] = ((-1//4 + 0//1*im)*eta_var^2) / ((1//1 + 0//1*im))
B[17,9] = ((0//1 + 1//4*im)*eta_var) / ((1//1 + 0//1*im))
B[18,9] = ((0//1 - 1//4*im)*eta_var) / ((1//1 + 0//1*im))
B[19,10] = ((-1//8 + 0//1*im)*eta_var^2) / ((1//1 + 0//1*im))
B[21,10] = ((0//1 - 1//4*im)*eta_var) / ((1//1 + 0//1*im))
YB[12] = ((1//1 + 0//1*im)) / ((1//1 + 0//1*im))
```
(mind to use eta and not eta\_var, while recopying, if you do. And add the \* )
B \ YB gives a bug :
```julia
ERROR: MethodError: no method matching abs2(::RationalPoly{Polynomial{…}, Polynomial{…}})
The function `abs2` exists, but no method is defined for this combination of argument types.
```
```julia
Stacktrace:
[1] _qreltype(::Type{RationalPoly{Polynomial{DynamicPolynomials.Commutative{…}, Graded{…}, Complex{…}}, Polynomial{DynamicPolynomials.Commutative{…}, Graded{…}, Complex{…}}}})
@ LinearAlgebra ~/.julia/juliaup/julia-1.11.2/share/julia/stdlib/v1.11/LinearAlgebra/src/qr.jl:341
[2] qr(A::Matrix{RationalPoly{Polynomial{…}, Polynomial{…}}}, arg::ColumnNorm; kwargs::@Kwargs{})
@ LinearAlgebra ~/.julia/juliaup/julia-1.11.2/share/julia/stdlib/v1.11/LinearAlgebra/src/qr.jl:424
[3] qr(A::Matrix{RationalPoly{Polynomial{…}, Polynomial{…}}}, arg::ColumnNorm)
@ LinearAlgebra ~/.julia/juliaup/julia-1.11.2/share/julia/stdlib/v1.11/LinearAlgebra/src/qr.jl:422
[4] \(A::Matrix{RationalPoly{Polynomial{…}, Polynomial{…}}}, B::Vector{RationalPoly{Polynomial{…}, Polynomial{…}}})
@ LinearAlgebra ~/.julia/juliaup/julia-1.11.2/share/julia/stdlib/v1.11/LinearAlgebra/src/generic.jl:1134
```
for the 0 elments, you should use :
```julia
0 + 0*eta
```
I hope it’s okay for you 🙂
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [January 3, 2025, 1:06pm UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/15 "2025-01-03T13:06:37Z")
> [@gustave\_Robichon](#):
>
> And add the \*
Could you edit your message to use code block formatting? Instead of this:
```md
> code
```
… do this:
````md
```julia
code
```
````
That will preserve the `*` symbols. EDIT: thanks!
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [January 3, 2025, 3:37pm UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/16 "2025-01-03T15:37:58Z")
> [@gustave\_Robichon](#):
>
> B \ YB gives a bug
I realize now this isn’t a bug, it’s just `\` behaving as documented. It’s doc string says that overdetermined (actually all non-square) systems are solved via:
> the minimum-norm least squares solution computed by a pivoted QR factorization
… so not what you want.
I suggest trying with Nemo.jl.
---
### Author:  [@thofma](https://discourse.julialang.org/u/thofma)
#### Post date: [January 4, 2025, 5:48pm UTC](https://discourse.julialang.org/t/how-to-exactly-solve-large-overdetermined-systems-of-linear-equations/86710/17 "2025-01-04T17:48:02Z")
To elaborate on what @nsjako wrote. Here is how you can do it using Nemo.jl:
```julia-repl
julia> QQi = Nemo.GaussianRationals();
julia> QQieta, eta_var = QQi[:eta];
julia> iim = QQi(0, 1);
julia> B = zero_matrix(QQieta, 21, 10);
julia> YB = zeros(QQieta, 21);
julia> begin
B[1,2] = ((-1//8 + 0//1*iim)*eta_var^2) / ((1//1 + 0//1*iim))
B[3,2] = ((0//1 + 1//4*iim)*eta_var) / ((1//1 + 0//1*iim))
B[4,3] = ((-1//4 + 0//1*iim)*eta_var^2) / ((1//1 + 0//1*iim))
B[5,3] = ((0//1 - 1//4*iim)*eta_var) / ((1//1 + 0//1*iim))
B[6,3] = ((0//1 + 1//4*iim)*eta_var) / ((1//1 + 0//1*iim))
B[7,4] = ((-1//8 + 0//1*iim)*eta_var^2) / ((1//1 + 0//1*iim))
B[8,4] = ((0//1 + 1//4*iim)*eta_var) / ((1//1 + 0//1*iim))
B[10,1] = ((1//1 + 0//1*iim)) / ((1//1 + 0//1*iim))
B[11,6] = ((-1//4 + 0//1*iim)) / ((1//1 + 0//1*iim))
B[12,7] = ((-1//4 + 0//1*iim)) / ((1//1 + 0//1*iim))
B[13,8] = ((-1//8 + 0//1*iim)*eta_var^2) / ((1//1 + 0//1*iim))
B[14,8] = ((0//1 - 1//4*iim)*eta_var) / ((1//1 + 0//1*iim))
B[16,9] = ((-1//4 + 0//1*iim)*eta_var^2) / ((1//1 + 0//1*iim))
B[17,9] = ((0//1 + 1//4*iim)*eta_var) / ((1//1 + 0//1*iim))
B[18,9] = ((0//1 - 1//4*iim)*eta_var) / ((1//1 + 0//1*iim))
B[19,10] = ((-1//8 + 0//1*iim)*eta_var^2) / ((1//1 + 0//1*iim))
B[21,10] = ((0//1 - 1//4*iim)*eta_var) / ((1//1 + 0//1*iim))
YB[12] = QQieta(((1//1 + 0//1*iim)) / ((1//1 + 0//1*iim)))
end;
julia> solve(B, YB; side = :right)
10-element Vector{AbstractAlgebra.Generic.Poly{Nemo.QQiFieldElem}}:
0
0
0
0
0
0
-4
0
0
0
```
If you want to go back and forth between `Complex{Rational{BigInt}}` and the Nemo type, you can do
```julia
julia> a = QQi(2, 3);
julia> Complex(Rational(real(a)), Rational(imag(a)))
2//1 + 3//1*im
```