Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

generic lufact requires more than it should #431

Closed
arghhhh opened this issue May 23, 2017 · 4 comments · Fixed by JuliaLang/julia#44571
Closed

generic lufact requires more than it should #431

arghhhh opened this issue May 23, 2017 · 4 comments · Fixed by JuliaLang/julia#44571

Comments

@arghhhh
Copy link

arghhhh commented May 23, 2017

I've made a minimal Mod{N} with only the operations required to be a field.
I needed to make a few modifications to the linalg code to get A \ b to work for a square matrix A with Mod{N} as the element type.

Specifically:
In a few places, a literal "0" is used instead of zero( .. ). I'd rather not add a conversion / promotion. It should work without this.
In one place "real( .. )" is used to initialize a variable that will hold the result of abs( .. )

If I turn off pivoting to reduce the requirements on my type further (so I don't need < or abs) then it doesn't work due to a SingularException error. I added a search for the first non-zero index and pivot on that - so its now pivoting, but without the requirement to have "abs" and "<".
I use this to line to turn off pivoting for matrices of my type Mod{N}:
import Base: lufact lufact{N}(A::AbstractMatrix{Mod{N}}, pivot::Union{Type{Val{false}}, Type{Val{true}}} = Val{false}) = lufact!(copy(A), pivot)

Below is a complete diff of my changes relative to v0.6.0-rc2.

Its pretty amazing that you can do this in Julia. It'd be great if we could fix these minor issues.
I'm not a mathematician, so I'm not completely sure about the pivoting on non-zero element, but it worked for me on my one(!) test case.


diff --git a/base/linalg/generic.jl b/base/linalg/generic.jl
index 6f44208..21e032e 100644
--- a/base/linalg/generic.jl
+++ b/base/linalg/generic.jl
@@ -955,7 +955,7 @@ true
 function istriu(A::AbstractMatrix)
     m, n = size(A)
     for j = 1:min(n,m-1), i = j+1:m
-        if A[i,j] != 0
+        if A[i,j] != zero( eltype(A) )
             return false
         end
     end
@@ -990,7 +990,7 @@ true
 function istril(A::AbstractMatrix)
     m, n = size(A)
     for j = 2:n, i = 1:min(j-1,m)
-        if A[i,j] != 0
+        if A[i,j] != zero( eltype(A) )
             return false
         end
     end
diff --git a/base/linalg/lu.jl b/base/linalg/lu.jl
index 7105bf3..ac2cf7f 100644
--- a/base/linalg/lu.jl
+++ b/base/linalg/lu.jl
@@ -39,7 +39,7 @@ function generic_lufact!(A::StridedMatrix{T}, ::Type{Val{Pivot}} = Val{true}) wh
             # find index max
             kp = k
             if Pivot
-                amax = real(zero(T))
+                amax = abs(zero(T))
                 for i = k:m
                     absi = abs(A[i,k])
                     if absi > amax
@@ -47,9 +47,16 @@ function generic_lufact!(A::StridedMatrix{T}, ::Type{Val{Pivot}} = Val{true}) wh
                         amax = absi
                     end
                 end
+        else # find first index with non-zero
+                for i = k:m
+                    if A[i,k] != zero(T)
+                        kp = i
+                        break
+                    end
+                end
             end
             ipiv[k] = kp
-            if A[kp,k] != 0
+            if A[kp,k] != zero(eltype(A))
                 if k != kp
                     # Interchange
                     for i = 1:n
@dpsanders
Copy link

Could you please submit a Pull Request with these changes so that people can easily test them out?
It would also be great to include some tests, for example with your ModN type.

@arghhhh
Copy link
Author

arghhhh commented May 30, 2017

OK. I am new to this (never contributed to any open source project).
There are really two issues here - one is relatively trivial uses of zero() and one() instead of 0 and 1 literals and the use of abs(zero(T)) instead of real(zero(T)) to avoid having to define real() which might not make sense for some matrix eltypes (eg for modulo{polynomial} ).

I put an example derived from the ModInt example in the Julia source here:
https://gist.github.com/arghhhh/abdb8b6d28039f683f5287e124903670
and will submit a pull request shortly.

The second issue is how to control pivoting in lufact. It seems to me that there could be three options here:

    • No pivoting
    • Non-zero pivoting - ensure that the pivot is non-zero, but don't attempt to choose largest
    • Pivoting - requires abs and < which might not make sense for some eltypes

The current implementation uses Val{false} and Val{true} to choose between 1. and 3. Adding 2. would require a bigger discussion.

@andreasnoack
Copy link
Member

It would be great with the other pivoting option as well. The reason for the current design is that we try to use the same API for LU and QR. The pivoted QR returns a different type which is the reason why we use Vals here. For LU, we use the same type so Val is not really needed but we should probably be consistent across LU and QR. We could consider adding pivoting types instead or use symbols instead of Booleans in the Val, e.g. Pivot{:anyrow}, Val{:maxrow}, or NormCol. I'm not super happy with any of these. Thoughts?

@arghhhh
Copy link
Author

arghhhh commented Jun 2, 2017

I had wondered why the Val{} argument had been used. I was thinking that using a symbol instead of the Bool would work - but I only realized this needed to be a separate option after my original issue submission. My original proposal hijacked Pivot=Val{false} to mean pivot on !iszero.

I'm not sure how to proceed. I'm new to contributing to Julia, so it was exciting(!) to see my PR make progress - but I haven't actually solved my original problem. My motivation is that I'm learning about finite fields and abstract algebra to enable understanding and implementation of error control coding - similar to https://github.com/andrewcooke/IntModN.jl

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

4 participants