Quantcast

Why doesn't lufact return Triangular Matrices for its data?

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Why doesn't lufact return Triangular Matrices for its data?

Chris Rackauckas
Why doesn't lufact return Triangular Matrices?

A = rand(1000,1000)
B = lufact(A,Val{false})
C = B[:U]
D = Base.LinAlg.UpperTriangular(B[:U])

The returned C is just a normal 1000x1000 array, not the same as D. I know there isn't a storage advantage to using Triangular matrices, but it seems there's a speed benefit to using Triangular matrices (for example I checked D*D instead of C*C), and maybe further improvements to Triangular matrices could be made after the array buffer implementation. 
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Why doesn't lufact return Triangular Matrices for its data?

Andreas Noack
The problem is that our Triangular types are square and U might not be square, e.g.

julia> lufact(randn(3,4))[:U]
3×4 Array{Float64,2}:
 -2.98058  -0.234937   1.49172   1.7675
  0.0      -1.15066   -1.72942   1.26475
  0.0       0.0       -1.58043  -0.0808058

We could loosen then type to allow for trapezoid matrices but it would make them more complicated and require more checks.

On Tue, Sep 27, 2016 at 12:08 PM, Chris Rackauckas <[hidden email]> wrote:
Why doesn't lufact return Triangular Matrices?

A = rand(1000,1000)
B = lufact(A,Val{false})
C = B[:U]
D = Base.LinAlg.UpperTriangular(B[:U])

The returned C is just a normal 1000x1000 array, not the same as D. I know there isn't a storage advantage to using Triangular matrices, but it seems there's a speed benefit to using Triangular matrices (for example I checked D*D instead of C*C), and maybe further improvements to Triangular matrices could be made after the array buffer implementation. 

Loading...