Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Suppose I have a $m\times n$ matrix $A$ (real for simplicity). Then SingularValueDecomposition[A] yields 3 matrices $U$, $\Sigma$ and $V$ such that

$A = U\Sigma V^\top = u_1 \sigma_1 v_1^\top + u_2 \sigma_2 v_2^\top + \cdots$,

where $U = [u_1\;\; u_2 \ldots]$ (vector = column matrix).

Successive approximations for $A$ are given by the rank 1 matrices $u_1 \sigma_1 v_1^\top$, etc. I wanted to make a function to compute this but found myself doing all sorts of messy manipulations due to Mathematica's matrix structure (and/or my lack of knowledge on how to use them).

Here is my tentative function. It has as input any matrix A and an optional argument n saying how many rank 1 matrices should be summed. This function is not meant to be numerically fast or anything; perhaps I would use it as a educational tool to examine different approximations of A or etc. The point is in the Sum function where I would like to know: is there a more efficient (or perhaps cleaner) way of doing the matrix multiplications?

RankOneApprox[A_, n_: 1] := Block[{U, Sig, V},

{U, Sig, V} = SingularValueDecomposition[A];

Sum[Sig[[i,i]] ({U[[All, i]]}\[Transpose].{V[[All, i]]}), {i, 1, n}]];
share|improve this question
add comment

1 Answer

up vote 8 down vote accepted

How about the following:

RankOneApprox[A_, n_: 1] := Module[
  {U, Sig, V},
  {U, Sig, V} = SingularValueDecomposition[A];
  U.DiagonalMatrix[Diagonal[Sig[[1 ;; n]]], 0, Length[a]].ConjugateTranspose[V]
  ]

Here I take the first n diagonal elements of Sig. That's equivalent to summing the projectors you wrote out explicitly. With these n singular values, I then form a new matrix that has zeros everywhere else, and insert that into the definition of the SVD. When n equals the dimension of a, you get the original SVD back.

This assumes a to be a square matrix - if desired, it would be pretty easy to extend to rectangular matrices. Here I just want to illustrate how one uses the matrix manipulations to avoid writing explicit sums. Also, your original function has the Transpose in a syntactically incorrect spot, and I've fixed that in the above version.

Edit

As mentioned by R.M in the comment, there is another way of calling SingularValueDecomposition that already weeds out the smallest singular values for you: SingularValueDecomposition[A, n]

However, this returns condensed versions of all the matrices, see also the documentation, where an example is given (under Applications) that does something very similar to this question. One has to be a little careful about the dimensionality of the factors in the SVD formula, because the unitary matrices in the condensed form are no longer represented by square matrices (in particular, the syntax Dot @@ SingularValueDecomposition in the comment won't work because conjugate-transposition was omitted). Fortunately, to construct the corresponding matrix approximation, we can maintain the same syntax as in the above function:

RankOneApprox[A_, n_: 1] := Module[
  {U, Sig, V},
  {U, Sig, V} = SingularValueDecomposition[A, n];
  U.Sig.ConjugateTranspose[V]
  ]
share|improve this answer
1  
It might be simpler to use the second argument of SVD to get the n largest singular values. Oneliner: Dot @@ SingularValueDecomposition[A, n] –  rm -rf Jul 14 '12 at 23:59
    
@R.M Nice simplification! –  Jens Jul 15 '12 at 0:24
    
@R.M OK- I added a version of your comment, with the addition of a missing transpose. –  Jens Jul 15 '12 at 0:56
    
Ah, yes of course :) My test matrix happened to have symmetric (and identical) u and v, so that was hasty. Nevertheless, thanks :) –  rm -rf Jul 15 '12 at 1:01
    
This looks great. Thank you very very much for the replies. –  Gabriel Landi Jul 15 '12 at 3:01
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.