Linear-Equation-System Solver("\") does not work reliably for the "Rational" Type

[Kurolong]

Hey Guys and Gals  ^^
I'm having trouble with Julia. Maybe one of you can help me out
The program i am writing requires the linear-equation-system-solver(command "\") in the folder "Julia-0.4.5\share\julia\base\linalg".
Now my problem is that i need the equation solved strictly on integers, an approximate solution is of no use to me, which is why i tried working with the "Rational" Type.

Curiously, the solver seems to work for the Rational Type if and only if no pivoting is required. Otherwise, the following Error is thrown:

WARNING: pivoting only implemented for Float32, Float64, Complex64 and Complex128
ERROR: LoadError: OverflowError()
 in * at rational.jl:188
 [inlined code] from linalg/generic.jl:471
 in reflector! at no file:0
 in A_ldiv_B! at linalg/qr.jl:346
 in \ at linalg/qr.jl:398
 in \ at linalg/dense.jl:450
 in include at boot.jl:261
 in include_from_node1 at loading.jl:320
while loading C:\users\<username omitted>\documents\julia-0.4.5\sg_project\v0.3\Test02.jl, in expression starting on line 3

Code to reproduce the Error:

A=[[2//1,0//1] [3//1,3//2] [1//1,3//2]]
v=[2//1,0//1]
A\v

Now i am under the impression that i could build a relatively simple workaround if i knew exactly where in the developer's code the actual problem is, but i am confused by the whole organization of the thing. Building my own solver could certainly be done, but this will most likely be ineffecient. Unfortunately, runtime is very important here, as the method i am currently devising is supposed to be embedded in some very expensive loops.

I am a very unexperienced programmer and will likely not understand much programmertalk, but i am a mathematician and will probably have little problems in this respect.

Maybe you have some ideas how to handle this problem.

Thank you in advance for taking the time to look at my problem :)

[Mauro]

This should open the editor in the right place:

julia> @edit A\v

If that does not work do this to look at it:

julia> @less A\v

Best probably monkey-patch it to work out the details, this is porgammer
speak for overloading the called method Base.(\)(A::StridedMatrix, B::StridedVecOrMat).

Once it works insert it back into dense.jl, rebuild, add tests, make a
pull request, voila!

[Chris Rackauckas]

You're not going to get good runtimes with this if you're using Rationals. Each time a command is done, it needs to reduce the fraction. If you're working with small examples like you show here, then yes this can be worthwhile. However, even with good programming I don't think would scale to something like 1000x1000 systems.

That said, this OverflowError() is may be because you're overflowing the integers? Try switching to Rational{BigInt} and see if you still get this error. This will decrease the performance even more though.

I know that being Mathematician the prospect of getting an exact answer sounds nice, but there's a reason why it's not done very often. It's can become orders of magnitude more computationally intensive than even precise approximations. If you need something which is really precise, try BigFloats or Decimals.jl.

[Kurolong]

The thing is, i did not write that i need an exact integer solution because i don't like approximation, but because i am working in N^n and rounding can cause catastrophically false results. The matrices i am working on are, however, not all that big (i would guess absolute max 30x100, usually more like 5x20).

[Kurolong]

Unfortunately, switching to Rational{BigInt} did not help. same error :(

[Kurolong]

Turns out that it actually did fix the error, but i was too tired and confused things.
BUT the components of the fractions were a lot too large.
A Gauß-algorithm should have returned very managable fractions, not fractions of 300-decimal-place-numbers by 300-decimal-place-numbers.
I suspect this is somehow caused by the lack of pivoting.
Maybe this somehow results in dividing by zero or extremely small values?

[Dan]

You can also try to use the pseudo-inverse $A^{+}$ in LaTeX. It is also called the Moore-Penrose inverse.

In this case (with variable names in thread and `AP` used for the inverse) it is:

 

AP = A'*inv(A*A')

And then we have ('x' is the solution to the equations):

x = AP*v

and as it should be

A*x == v

is true.

Hope this helps.

-- Dan

[Kurolong]

But Julia cannot calculate the Moore-Penrose-Inverse for the Rational Type :(

[Steven G. Johnson]

On Friday, July 15, 2016 at 10:31:20 PM UTC+1, Kurolong wrote:

Turns out that it actually did fix the error, but i was too tired and confused things.
BUT the components of the fractions were a lot too large.
A Gauß-algorithm should have returned very managable fractions, not fractions of 300-decimal-place-numbers by 300-decimal-place-numbers.

How do you know that?   In general, the number of digits in the denominators for exact rational arithmetic increases linearly with the number of computational steps.   Since each output in A \ b by Gaussian elimination on an nxn matrix A requires O(n^2) arithmetic operations (O(n^3) operations for n outputs), 300 decimal places is reasonable even for a 10x10 matrix.

(This is why rational bignum arithmetic is usually even more of an impractical luxury than ordinary bignum arithmetic.)

[Kurolong]

Well, firstly, with matrices, the increase should be distributed among entries.
secondly, i was talking about running my above example in BigInt.
A 2x3 matrix should most definitely not produce fractions like that.
Thirdly, the system has the solution [0.8,-2.0,0.4].
This is not what the function returned.

[Andreas Noack]

A couple of things go wrong here. Right now, Julia tries to use the QR factorization to solve underdetermined systems. If things had worked the way I'd planned, your matrix would have been promoted to floating point elements (I know that it's not what you want so keep reading). The promotion fails and Julia tries to QR factorize your matrix even though it has rational elements. What happens afterwards is still a bit unclear to me but I'm looking into it.

We should add an LU based solver for underdetermined systems which should get called for rationals. Unfortunately due to our release cycle, it might take a while to get into a release version of Julia. Meanwhile something like

julia> A = Matrix{Rational{BigInt}}(rand(1:10, 4, 5))

4x5 Array{Rational{BigInt},2}:

 6//1  9//1   6//1  5//1  3//1

 6//1  2//1  10//1  9//1  9//1

 4//1  1//1   1//1  4//1  1//1

 6//1  3//1   6//1  7//1  5//1

julia> b = A*ones(Int, size(A, 2));

julia> F = lufact(A);

julia> x = copy!(zeros(eltype(b), size(A, 2)), b);

julia> A_ldiv_B!(UpperTriangular(sub(F[:U], :, 1:size(A, 1))), A_ldiv_B!(LowerTriangular(F[:L]), sub(x, 1:size(A, 1))));

julia> A*x - b

4-element Array{Rational{BigInt},1}:

 0//1

 0//1

 0//1

 0//1

shouldn't be too inefficient (if wrapped in a function).

On Friday, July 15, 2016 at 2:50:01 PM UTC-4, Kurolong wrote:

[Kurolong]

Thank you, this was very helpful.
Just a small question: Are the UpperTriangular and LowerTriangular commands actually necessary in some case?
The other question i have is what A_div_B! actually does in this case, but i have a feeling that this is hard to explain.

[Andreas Noack]

The UpperTriangular and LinAlg.UnitLowerTriangular parts are not necessary. Julia will detect that but since we already know that they are triangular we can save the costs of the checks.

A_ldiv_B!(A,b) is indeed ugly. It's just A\b where the result is stored in b. That is useful here because it will save us a couple of temporary allocations.

[Kurolong]

Ah right, so what you are doing is just solving successively for L and U.
Okay, all questions answered.
Thank you^^