Posts Tagged ‘matrix’

« Older Entries

Setting row of sparse matrix to zero in matlab

Friday, May 24th, 2013

I wanted to set a whole row of a sparse matrix A to zero:


A = sparse(1e8,1e8);
tic;
A(r,:) = 0;
toc

But this was rather slow. About 1.7 seconds.

Turns out all rows and columns are not created equal in matlab. Because of how matlab stores sparse matrices it’s much faster to set a whole column to zero:


tic;
A(:,c) = 0;
toc

This is about 0.0002 seconds.

So if you can refactor your code you’re better off messing with columns. If you can’t refactor, you might think it would be fast to transpose before and after setting a row to zero:


tic;
A = A';
A(:,r) = 0;
A = A';
toc

But this helps not at all. It seems transposes are just stored as flags, which is usually a good thing. But here it means we don’t gain anything, since after the transpose, now rows the cheap route. This also implies that refactoring could be nontrivial.

One might also think you could first find all the non-zero entries and then set them each to zero. But just the rowwise find is expensive:


tic;find(A(10,:));toc

~1 second

In fact just the row-wise access is expensive


tic;A(10,:);toc

~1 second

Seems like there’s no happy answer.

Tags: , , , , ,
Posted in code | No Comments »

Initializing Eigen matrix to zeros, the wrong way.

Tuesday, April 30th, 2013

I recently found a pretty silly bug in the way I was initializing an Eigen matrix to be all zeros. Say, my matrix was A, then I wrote this (admittedly in haste) to resize to an m by n matrix of zeros using:


A.resize(m,n);
A *= 0;

This worked most of the time. As should be expected if A is initialized to numbers. But unfortunately, the resize doesn’t initialize the values to anything, and some of the random leftovers in memory will be, when interpreted as floats, not-a-numbers (NaNs). Thus when I thought I was zeroing them out I was keeping the NaNs: 0*NaN = NaN. Luckily, this is also why NaN-related bugs are so easy to trace.

Of course the correct thing to do is explicit set each value to zero. Eigen even has a dedicated function:


A.setZero(m,n);

Tags: , , , , ,
Posted in code | No Comments »

Fast, sparse kronecker product with identity in MATLAB

Friday, July 6th, 2012

I needed to compute:


  K = kron(A,speye(n,n));

but was a little dismayed at the speed. Here’s how I get about a 6x speed up:


  % Compute kron(speye(n,n),A) in a fast way       
  C = cell(n,1);                                   
  [C{:}] = deal(A);                                
  K = blkdiag(C{:});                               
  % Rearrange rows and columns                     
  I = reshape(reshape(1:(size(A,1)*n),size(A,1),n)',size(A,1)*n,1);
  J = reshape(reshape(1:(size(A,2)*n),size(A,2),n)',size(A,2)*n,1);
  K = K(I,J);

Tags: , , , , , , ,
Posted in code | No Comments »

LU decomposition using two Cholesky decompositions (part two)

Sunday, March 4th, 2012

I previously posted about how often you can compute a LU decomposition of a matrix Q using two Cholesky decompositions. Specifically when your matrix is of the form:


Q = / A  C \
    \ CT * /

In that post I assumed you had a subroutine that computes the Cholesky factorization of a symmetric positive-definite matrix A:


LLT=A

Where L is a lower-triangular matrix.

Often, routines exist to give you an even better factorization which includes a permutation matrix S:


LLT=STAS

This often results in a more sparse L.

Working through the LU decomposition using this type of Cholesky factorization is a bit tedious and requires some careful book keeping.

This time I want to fill in the missing parts of the following system:


/ ? * \ / A  C \ / ? * \ = / ? * \ / ? ? \
\ * ? / \ CT * / \ * ? /   \ ? ? / \ * ? /

Since A is symmetric positive definite we use our factorization STAS = LLT which gives us:


/ ST * \ / A  C \ / S * \ = / L * \ / LT ? \
\ *  ? / \ CT * / \ * ? /   \ ? ? / \ *  ? /

But this means that know we need that STC = L ?. But L is lower triangular so it’s easy to compute:


M = L-1STC

We can then place M into our LU decomposition:


/ ST * \ / A  C \ / S * \ = / L  * \ / LT M \
\ *  ? / \ CT * / \ * ? /   \ MT ? / \ *  ? /

Now at this point we could use cholesky factorization without the permutation matrix to complete the missing parts, but we might as well use our better decomposition again. To do this we must be a bit careful. We want to create the zeros in the lower right corner of the system. To do this we need a factorization of MTM. But our new routine gives us:


KKT = TTMTMT

where K another lower triangular matrix and T another permutation matrix.

So we replace M with MT in our system and plug in K. This gives us:


/ ST * \ / A  C \ / S * \ = / L    * \ / LT  MT \
\ *  ? / \ CT * / \ * ? /   \ TTMT K / \ *  -KT /

Finally, compute the missing pieces by evaluating the right hand side and we indeed get our original system conjugated with a permutation matrix composed of the permutation matrices from both of the cholesky factorizations:


/ ST *  \ / A  C \ / S * \ = / L    * \ / LT  MT \
\ *  TT / \ CT * / \ * T /   \ TTMT K / \ *  -KT /

Tags: , , , , , , , , ,
Posted in code | No Comments »

2D rotation matrix plus another 2D rotation matrix is a similarity matrix (scaled 2D rotation)

Wednesday, November 16th, 2011

Consider we have two rotation matrices:


R1 = / cos(θ1) -sin(θ1) \
     \ sin(θ1)  cos(θ1) /

R2 = / cos(θ2) -sin(θ2) \
     \ sin(θ2)  cos(θ2) /

Then we can show that:


R1 + R2 = s * R3

where s is a scalar and R3 is another 2D rotation.

First we just add the components to get:


R1 + R2 = / cos(θ1)+cos(θ2)   -(sin(θ1)+sin(θ2)) \
          \ sin(θ1)+sin(θ2)    cos(θ1)+cos(θ2)   /

Notice we already have the right pattern:


R1 + R2 =  / A -B \
           \ B  A /

we just need to show that we can pull out a common scaling term and angle from A and B.

We can use trigonometric identities to turn the sums of cosines and the sums of sines into products, namely:


cos(θ1)+cos(θ2) = 2*cos((θ1+θ2)/2)*cos((θ1-θ2)/2)
sin(θ1)+sin(θ2) = 2*sin((θ1+θ2)/2)*cos((θ1-θ2)/2)

Now we declare that:


s = 2*cos((θ1-θ2)/2)
θ3 = (θ1+θ2)/2

Then it follows that:


R1 + R2
  =
  / s*cos((θ1+θ2)/2)  -s*sin((θ1+θ2)/2) \
  \ s*sin((θ1+θ2)/2)   s*cos((θ1+θ2)/2) /
  =
  s * / cos(θ3)  -sin(θ3) \
      \ sin(θ3)   cos(θ3) /
  = 
  s * R3

Tags: , , , , ,
Posted in code | No Comments »

LU decomposition using two Cholesky factorizations (part one)

Thursday, November 3rd, 2011

My colleague showed me a neat linear algebra trick to compute the LU decomposition of a special type of matrix using two Cholesky factorizations. This is typical much faster.

We’d like to factorize Q = LU where L is lower triangular and U is upper triangular. If Q is symmetric positive-definite we can use Cholesky factorization to compute Q = LLT. Often this unfortunately is not the case and we typically fall back on the slower LU decomposition. But we can use a trick if Q is a special matrix of the following block form:


Q = / A  C \
    \ CT * /

where A is an n by n symmetric, positive definite matrix, C is a full row rank m by n matrix and * is all zeros (here m by m).

These types of systems are actually very common. They arrise often during energy minimization with linear equality constraints. For example I might have the quadratic energy:


minimize xT A x + xT B
subject to CT x= D

Using lagrange multiplers I can transform this minimization into finding the sandle point of


[xT λT] / A  C \ /x\ + [xT λT] /  B  \
        \ CT * / \λ/           \ -2D /

And to find the sandle point we take derivatives with respect to each x and λ resulting in 


/ A  C \ /x\ = / -B/2 \
\ CT * / \λ/   \  D   /

Then we can solve for x and λ by computing an LU factorization of 


/ A  C \
\ CT * /

Just like our matrix is divided into 4 blocks, we can divide our target L and U into 4 blocks:


/ A  C \ = / ?  * \ / ?  ? \
\ CT * /   \ ?  ? / \ *  ? /

Since we are assuming that A is symmetric positive definite we can compute a Cholesky factorization, namely let A = LLT. Then we can place L into the appropriate corners of the lower and triangular parts:


/ A  C \ = / L  * \ / LT  ? \
\ CT * /   \ ?  ? / \ *   ? /

But this means that know we need that C = L ?. But L is lower triangular so it’s easy to compute:

 
M = L-1C

We can then place M into our LU decomposition:


/ A  C \ = / L  * \ / LT  M \
\ CT * /   \ MT ? / \ *  ? /

We’re almost done. All that’s left is reproduce the zeros in the bottom right corner on the left hand side. This means we need the two missing pieces when multiplied to equal -MTM. But again this is made easy via Cholesky facotrization. MTM is clearly symmetric so we can Choleksy factorize it into:


MTM = KKT

where K is lower triangular.

So filling in with K we have our finally LU decomposition:


/ A  C \ = / L   * \ / LT  M \
\ CT * /   \ MT -K / \ *  KT /

We only needed two cholesky facotrizaions. And further each Cholesky factorization can be completed according to the characteristics of A (our original quadratic coefficients) and C our constraints. So if A is sparse then we can use sparse cholesky factorization to compute L. But then say if C on the other hand is dense, M will come out dense and we can use dense Cholesky factorization on MTM to compute K.

Update: See my following post about taken advantage of cholesky decomposition which give you permutation matrices.

Tags: , , , , , ,
Posted in code | No Comments »

Find minimum (or maximum, etc.) non-zero entry in a sparse matrix in MATLAB

Wednesday, November 2nd, 2011

Actually this should also work for dense matrices, that is it will give you the smallest entry in the matrix not equal to zero.


  min(A(find(A)))

Tags: , , ,
Posted in code | 1 Comment »

Eigen gotcha, sparse matrix inneriterator type

Sunday, October 9th, 2011

I just experienced a really annoying Eigen gotcha. If I have a sparsematrix A, defined as:


SparseMatrix<double> A;
// ... filled in somehow

then I can loop over the values with:


for(int k=0; k<A.outerSize(); ++k)
{
  // Iterate over inside
  for(SparseMatrix<double>::InnerIterator it (A,k); it; ++it)
  {
      // it.row(),  it.col(), it.value()
  }
}

This works fine. But if instead I use the wrong type for the InnerIterator, say int, I get no complaint from the compiler, just occasional garbage when I use it.col(), it.row(), or it.value():


for(int k=0; k<A.outerSize(); ++k)
{
  // Iterate over inside
  for(SparseMatrix<int>::InnerIterator it (A,k); it; ++it)
  {
      // it.row(),  it.col(), it.value() produce occasional garbage
  }
}

Tags: , , , , ,
Posted in code | No Comments »

Array multiply each column of sparse matrix by a column vector

Tuesday, August 16th, 2011

I was using essentially a dense matrix operation to multiply each column in a sparse matrix by a constant column vector:


S = sprandn(1000,1000,0.01);
% random sparse matrix
S = sprandn(1000,1000,0.01);
% random sparse matrix
S = sprandn(1000,5000,0.01);
% random column vector of same size
v = rand(size(S,1),1);
% old way: repeat column vector for each column in sparse matrix
old = S .* repmat(v,1,size(S,2));

It’s much faster to use the builtin bsxfun function with the binary operation times:


% new way: use bsxfun with binary times function
new = bsxfun(@times,S,v);

Update: There is a much better way, by utilizing sparse matrix-matrix multiplication:


new = diag(sparse(v))*S;

Source

Tags: , , , , ,
Posted in code | 2 Comments »

Common mesh energies, sum notation and matrix notation

Monday, June 20th, 2011

I’m always re-deriving the matrix notation of common triangle mesh energies like the Dirichlet energy, Laplacian energy, or local rigidity energy. Here’s my usual strategy. I’ll assume that if your energy came from a continuous integral you can at least get it to a sum over mesh elements. Here I write my sums over vertices and their neighbors, but surely you could also sum over triangles.

mesh energies sum to matrix form
Pretty-printed pdf document

Uglier plain text, unicode version:


∑i∑j∈N(i) vi  = ∑i * #N(i) * v = VT * D * 1, 
  where D is diagonal with Dii = #N(i)
∑i∑j∈N(i) vj  = VT * A * 1, where A is the adjacency matrix
∑i∑j∈N(i) vi - vj = VT * D * 1 - VT * A * 1 = VT * (D-A) * 1

∑i∑j∈N(i) wij  = 1 * Aw * 1, where Aw adjacency with Awij = wij, if j∈N(i)
∑i∑j∈N(i) wij*vi  = ∑i vi * ∑j∈N(i) wij  = VT * Dw * 1, 
  where Dw diagonal with Dwii = ∑j∈N(i) wij
∑i∑j∈N(i) wij*vj  = VT * Aw * 1

∑i mi * ∑j∈N(i) wij vi = ∑i mi * vi * ∑j∈N(i) wij = VT * M * Dw * 1,
  where M is diagonal with Mii = wi
∑i mi * ∑j∈N(i) wij vj = VT * Aw * M * 1

∑i∑j∈N(i) vi² = VT * D * V
∑i∑j∈N(i) vivj = VT * A * V
∑i∑j∈N(i) vj² = VT * RD * V, where RD is diagonal, RDjj = ∑i 1 if j∈N(i)
  if j∈N(i) implies i∈N(j) then
  D = RD
  so ∑i∑j∈N(i) vj² = ∑i∑j∈N(i) vi² = VT * D * V
∑i∑j∈N(i) (vi - vj)² = ∑i∑j∈N(i) vi² - 2 * ∑i∑j∈N(i) vivj + ∑i∑j∈N(i) vj²
                     = 2 * ∑i∑j∈N(i) vi² - 2 * ∑i∑j∈N(i) vivj
                     = 2 * VT * D * V - 2 * VT * A * V
                     = 2 * VT * (D - A) * V
                     = 2 * VT * L * V

∑i∑j∈N(i) wij*vi²  = ∑i vi² * ∑j∈N(i) wij = VT * Dw * V
∑i∑j∈N(i) wij*vivj = VT * Aw * V
∑i∑j∈N(i) wij*vj² = VT * RDw * V, where DW is diagonal, RDwjj = ∑i wij if j∈N(i)
  if j∈N(i) implies i∈N(j) and wij = wji then 
  Dw = RDw
  so ∑i∑j∈N(i) wij*vj² = ∑i∑j∈N(i) wij*vi² = VT * Dw * V
∑i∑j∈N(i) wij(vi - vj)² = 
  = ∑i∑j∈N(i) wij*vi² - 2 * ∑i∑j∈N(i) wij*vivj + ∑i∑j∈N(i) wij*vj²
  = 2 * ∑i∑j∈N(i) wij*vi² - 2 * ∑i∑j∈N(i) wij*vivj
  = 2 * VT * Dw * V - 2 * VT * Aw * V
  = 2 * VT * (Dw - Aw) * V
  = 2 * VT * Lw * V

∑i mi * ∑j∈N(i) wij vi² = ∑i mi * vi² * ∑j∈N(i) wij = VT * M * Dw * V
∑i mi * ∑j∈N(i) wij vj² = VT * RDw * M * V
∑i mi * ∑j∈N(i) wij vivj = VT * Aw * M * V
∑i mi * ∑j∈N(i) wij(vi - vj)² = 
  = ∑i mi * ∑j∈N(i) wij*vi² 
    - 2 * ∑i mi * ∑j∈N(i) wij*vivj 
    + ∑i mi *∑j∈N(i) wij*vj²
  = 2 * ∑i mi * ∑j∈N(i) wij*vi² - 2 * ∑i mi * ∑j∈N(i) wij*vivj
  = 2 * VT * M * Dw * V - 2 * VT * Aw * M * V
  = 2 * VT * M * (Dw - Aw) * V
  = 2 * VT * M * Lw * V
  Note: let wij = cwij = cot(αij) + cot(βij) / mi and let mi = voronoi area at 
  vi, then the familar cotangent form of the discrete Dirichlet energy is
  EDirichlet = 2 * VT * M * Lcw * V

Tags: , , , , ,
Posted in code | No Comments »

« Older Entries