JuliaLang / LinearAlgebra.jl

Julia Linear Algebra standard library
Other
17 stars 4 forks source link

lufact should allow more pivoting options? #205

Closed gwhowell closed 4 years ago

gwhowell commented 9 years ago

I don't see the option of column pivoting in lufact

available in UMFPACK and in the matlab and octave sparse lu

Using column pivoting often gives a sparser and better conditioned L than just using row pivoting.

tkelman commented 9 years ago

I agree this would be quite useful. Another entry to my list in https://github.com/JuliaLang/julia/pull/10402#issuecomment-77775800 - pivoting is a more complicated subject than a boolean true/false. cc @andreasnoack

andreasnoack commented 9 years ago

I've spent most the my time in SuiteSparse land inside CHOLMOD and SPQR so I'm not familiar with all the options in UMFPACK. I've just browsed through the sources and also looked at the MATLAB documentation and from that it is not clear to my what is requested here. @gwhowell Could you please point to a place where this is described. Preferely in the UMFPACK source or manual, but just explaining how you'd do this in MATLAB would be helpful.

gwhowell commented 9 years ago

Hi,

My current interest is in rectangular A, then instead of the normal

equations

A^t A x = A^t b

I can work with an equivalent equation

L^t L y = z

Since L is usually quite a bit better conditioned than A, convergence of conjugate gradient on L^t L tends to be quite a bit faster. I think that's an easy change to make in the code (as easy as just taking out the check for square matrices ? )

For sparse matrices the column (as well as row) pivoting L is often better conditioned.

Here's the octave "help lu", this syntax may be a bit cleaner than the matlab .. will look for that next, and perhaps code snippets that use them both.

It's great you're looking at this so quickly, wow !!

'lu' is a function from the file /usr/lib64/octave/3.6.4/oct/x86_64-redhat-linux-gnu/lu.oct

-- Loadable Function: [L, U] = lu (A) -- Loadable Function: [L, U, P] = lu (A) -- Loadable Function: [L, U, P, Q] = lu (S) -- Loadable Function: [L, U, P, Q, R] = lu (S) -- Loadable Function: [...] = lu (S, THRES) -- Loadable Function: Y = lu (...) -- Loadable Function: [...] = lu (..., 'vector') Compute the LU decomposition of A. If A is full subroutines from LAPACK are used and if A is sparse then UMFPACK is used. The result is returned in a permuted form, according to the optional return value P. For example, given the matrix 'a = [1, 2; 3, 4]',

      [l, u, p] = lu (A)

 returns

      l =

       1.00000  0.00000
       0.33333  1.00000

      u =

       3.00000  4.00000
       0.00000  0.66667

      p =

       0  1
       1  0

 The matrix is not required to be square.

When called with two or three output arguments and a spare input matrix, 'lu' does not attempt to perform sparsity preserving column permutations. Called with a fourth output argument, the sparsity preserving column transformation Q is returned, such that 'P * A Q = L \ U'.

 Called with a fifth output argument and a sparse input matrix, 'lu'
 attempts to use a scaling factor R on the input matrix such that 'P
 * (R \ A) * Q = L * U'.  This typically leads to a sparser and more
 stable factorization.

 An additional input argument THRES, that defines the pivoting
 threshold can be given.  THRES can be a scalar, in which case it
 defines the UMFPACK pivoting tolerance for both symmetric and
 unsymmetric cases.  If THRES is a 2-element vector, then the first
 element defines the pivoting tolerance for the unsymmetric UMFPACK
 pivoting strategy and the second for the symmetric strategy.  By
 default, the values defined by 'spparms' are used ([0.1, 0.001]).

 Given the string argument 'vector', 'lu' returns the values of P
 and Q as vector values, such that for full matrix, 'A (P,:) = L *
 U', and 'R(P,:) * A (:, Q) = L * U'.

 With two output arguments, returns the permuted forms of the upper
 and lower triangular matrices, such that 'A = L * U'.  With one
 output argument Y, then the matrix returned by the LAPACK routines
 is returned.  If the input matrix is sparse then the matrix L is
 embedded into U to give a return value similar to the full case.
 For both full and sparse matrices, 'lu' loses the permutation
 information.

On Fri, Apr 24, 2015 at 8:04 AM, Andreas Noack notifications@github.com wrote:

I've spent most the my time in SuiteSparse land inside CHOLMOD and SPQR so I'm not familiar with all the options in UMFPACK. I've just browsed through the sources and also looked at the MATLAB documentation and from that it is not clear to my what is requested here. @gwhowell https://github.com/gwhowell Could you please point to a place where this is described. Preferely in the UMFPACK source or manual, but just explaining how you'd do this in MATLAB would be helpful.

— Reply to this email directly or view it on GitHub https://github.com/JuliaLang/LinearAlgebra.jl/issues/205.

gwhowell commented 9 years ago

Hi Andreas,

Here's the matlab help lu, Actually the syntax looks identical to the octave .. In addtion to a column pivot vector q, they can also have a scaling matrix R -- specified by how many arguments you put on the left hand side of the call (an additional argument on the right hand side is whether to have threshold pivoting, i.e., only pivot if the threshold is violated). Will look next at the documentation for SuiteSparse.

http://www.mathworks.com/help/matlab/ref/lu.html

[edit: removed MathWorks documentation text, replaced it with a link – @StefanKarpinski]

andreasnoack commented 9 years ago

It is my impression (from the UMFPACK manual and MATLAB manual) that UMFPACK is defaulting to a combination of "symmetric" pivoting where both rows and columns are interchanged and traditional partial pivoting where the rows are interchanged.

To be specific, which input argument in MATLAB or Octave is activating the row pivoting that you are asking for?

jiahao commented 9 years ago

@gwhowell please try to respect The MathWorks's copyright to the Matlab documentation. It is unclear what exactly you are referring to in the enormous block of text you have just copied.

gwhowell commented 9 years ago

Hi, You already have row pivoting, what I need most is the rectangular matrices (should be an easy fix, just don't throw out the rectangular matrices with an if statement)

matlab and octave specify the row and column pivoting options by asking for them in the output list

[l,u,p] = lu(a)

returns l, u, p -- p row pivot matrix so that lu = p \ a

[l,u,p] - = lu(a,"vector")

gives a permutation vector

[l,u,p,q] returns matrices p and q such that

l * u = p * a * q

Have a meeting, will look at the UMFPACK interface, hopefully Monday.

Regards, Gary Howell

On Fri, Apr 24, 2015 at 1:46 PM, Andreas Noack notifications@github.com wrote:

It is my impression (from the UMFPACK manual and MATLAB manual) that UMFPACK is defaulting to a combination of "symmetric" pivoting where both rows and columns are interchanged and traditional partial pivoting where the rows are interchanged.

To be specific, which input argument in MATLAB or Octave is activating the row pivoting that you are asking for?

— Reply to this email directly or view it on GitHub https://github.com/JuliaLang/LinearAlgebra.jl/issues/205.

andreasnoack commented 9 years ago

I've changed the title of this issue. It would probably be good to support more pivoting options at some point, but it is not a top priority right now.

tkelman commented 9 years ago

I think this is still a good long-term feature request to keep open, any specific reason to close it already?

stevengj commented 4 years ago

Closing as a duplicate of a more modern issue: #36720