LAPACK Archives

### [Lapack] psuedo inverse algorithm

 ``` Hi Adam, Few comments. - Please see my post on the LAPACK forum: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=5&t=160 on how to compute the pseudo inverse with LAPACK. Essentially you can do this through the (1) QR factorization of the matrix A (DGEQRF), (2) you can also use partial pivoting (DGEQP3), (3) you can also use SVD. Of course (4) you can always go straight at the formula: ( A^T A )^(-1) A^T, but that will not be numerically stable if A is ill-conditioned. - LAPACK does not have a routine to compute the pseudo-inverse of a matrix. From time to time we wonder if we should have such a routine, No user has ever asked for a routine to compute the pseudo-inverse, (besides for the reason: "I want to compute the pseudo inverse, because I want to"), and I am not aware of any application which requires it. - The main advantage of (2, QR with pivoting) and (3, SVD) is that these methods enable filtering. So that you can discard small singular values of A before inverting. This feature is needed for most applications. - If you go with a QR based method (1 or 2), you number of flops would be: (assuming m>=n) QR factorization: 2 mn^2 - 2/3n^3 Construct Q: 2 mn^2 - 2/3n^3 R \ Q^T: mn^2 So cost is 5 mn^2 - 4/3 n^3. This counts multiplications and additions. You may divide by two for multiplications only. - Now if we go with ( A^T A )^(-1) A^T, (not stable but let us do this anyway,) A^TA : (SYRK ): mn^2 Cholesky factorization of A^T A : ( POTRF ): n^3 / 3 Cholesky solve ( A^T A )^(-1) A^T: 2mn^2 So cost is 3 mn^2 + 1/3 n^3. This counts multiplications and additions. You may divide by two for multiplications only. Now looking at the algorithm you are giving in your post: - The algorithm does not look stable. I think it is as stable as computing the pseudo-inverse through the formula ( A^T A )^(-1) A^T. Essentially if you compute || A^dagger * A - I || you will see a dependency in K( A )^2 where K( A ) is the condition number of the matrix A. A "stable" algorithm should get a decency in kappa(A). - The algorithm has less flops when using elementary transformation, I believe we end up in the same number of flops for a QR factorization based method when the algorithm uses Householder reflectors. LAPACK uses Householder reflectors for pretty much anything except LU factorization, Cholesky factorization, etc. For example, we could very well use elementary operations for reduction to Hessenberg form, and this would pretty much valid, and halves the number of flops. We use Householder reflectors for most of these operations (despite the double number of flops) because these transformations are always stable. So all in all, this is a quick review of the algorithm. It is fun. I was not aware about it. I am not sure we are interested in including it in LAPACK, the main reason is that (1) not sure there is a user base for such a functionality, (2) we can obtain the pseudo-inverse from functionality in LAPACK, (3) the proposed algorithm is not the most stable one. Cheers, Julien. On Mar 2, 2013, at 9:17 AM, adam W > wrote: Hello, I have recently published an algorithm that computes the Moore-Penrose pseudoinverse, and does so in exact arithmetic. It accomplishes this with 2mn^2 multiplications (dimensions with m>n), thus is comparable to a single matrix multiply. I believe your function makes use of singular value decomposition. Mine does not achieve the singular values, which is why it works quicker for just the pseudoinverse. Here is the (thus far only) link to my algorithm: http://math.stackexchange.com/a/317053/43193 I believe that this could enhance your wonderful software, how may I help this happen? Do I need to submit publication into a journal, or maybe write some code? (I have my code in python currently, but the algorithm is quite simple really and could easily be ported-- at least it would be easy for me to do.) -Adam _______________________________________________ Lapack mailing list Lapack@Domain.Removed http://lists.eecs.utk.edu/mailman/listinfo/lapack -------------- next part -------------- An HTML attachment was scrubbed... URL: http://lists.eecs.utk.edu/mailman/private/lapack/attachments/20130305/8470d0c1/attachment.html ```
 Current Thread [Lapack] psuedo inverse algorithm, adam W [Lapack] psuedo inverse algorithm, Langou, Julien <=

For additional information you may use the LAPACK/ScaLAPACK Forum.
Or one of the mailing lists, or