Inversion of matrix made of diagonal blocks

July 20th, 2014

Recently I needed to invert (a.k.a. solve against) a large matrix with a special form. It was a block matrix where each block was diagonal:

 M = [M11 M12 ... M1n
      M21 ...
      ...
      Mn1 Mn2 ... Mnn];

where each Mij was a diagonal matrix. In my case Mij=Mji, though I won’t rely on this for the following.

Note that this is not the same as a “block diagonal matrix” though it can be rearrange into such a matrix.

The inverse of my kind of matrix will have the same sparsity pattern. I can find it using a recursive application of the blockwise matrix inversion formula.

Here’s a little function to do exactly that:

function N = diag_blk_inv(M,k)
  % DIAG_BLK_INV Invert a block matrix made out of square diagonal blocks. M
  % should be of the form:
  %
  %   M = [M11 M12 ... M1n
  %        M21 ...
  %        ...
  %        Mn1 Mn2 ... Mnn];
  % where each Mij is a diagonal matrix.
  %
  % Inputs:
  %   M  n*k by n*k matrix made out of diagonal blocks.
  % Outputs:
  %   N  n*k by n*k matrix inverse of M
  %

  switch k
  case 1
    assert(isdiag(M));
    N = diag(sparse(1./diag(M)));
    return;
  end

  assert(size(M,1)==size(M,2),'Must be square');
  assert(rem(size(M,1),k)==0,'Must be divisble by k');
  n = size(M,1)/k;
  % Extract 4 blocks
  A = M(0*n+(1:(k-1)*n),0*n+(1:(k-1)*n));
  B = M(0*n+(1:(k-1)*n),(k-1)*n+(1:n));
  C = M((k-1)*n+(1:n),0*n+(1:(k-1)*n));
  D = M((k-1)*n+(1:n),(k-1)*n+(1:n));
  assert(isdiag(D));
  % https://en.wikipedia.org/wiki/Invertible_matrix#Blockwise_inversion
  Ainv = diag_blk_inv(A,k-1);
  % Schur complement 
  S = (D-C*Ainv*B);
  assert(isdiag(S));
  Sinv = diag_blk_inv(S,1);
  N = [Ainv + Ainv*B*Sinv*C*Ainv -Ainv*B*Sinv; ...
       -Sinv*C*Ainv              Sinv];
end

Update: No doubt someone will see that I’m computing the inverse explicitly and call me a heretic. Actually in the case of diagonal matrices there’s some room for argument why one might want to do this. In terms of efficiency, this is orders of magnitude faster than calling inv(M)*B or M\B. However, matlab’s \ doesn’t seem to figure out that it could use cholesky decomposition. In that case L = chol(M);L\(L'\B) is roughly the same speed as diag_blk_inv(M,k)*B: both are fast.

Update: To be my own policeman, it probably is safer to use chol in this case. Otherwise I would in general I would need to check the scaling of M. It’s unfortunate though because a priori I don’t have any reason to believe chol will see the special structure of this matrix.

Interleave rows of single matrix

July 18th, 2014

Dealing with xyz coordinates in linear algebra for computer graphics, there are predominant two styles when concatenating a list of many positions/vectors/matrices. The first groups common coordinates together. For example a list of mesh vertex coordinates may look like:

V = [x1
     x2
     ...
     xk
     y1
     y2
     ...
     yk
     z1
     z2
     ...
     zk];

The other style groups common objects together, listing all coordinates in a stream. For our vertex positions this would look like:

V = [x1
     y1
     z1
     x2
     y2
     z2
     ...
     xk
     yk
     zk];

In OpenGL and array language the first option has a data access stride of 1 if pulling x-coordinates and the second option has a stride of n. Of course if pulling position vectors then the first option has a stride of n and the second a stride of 1. Hence, the two styles really have appropriate use cases depending on data access patterns.

This gets a bit more complicated if we’re not just storing a scalar in each row, but a n-vector.

In matlab, it’s very easy to convert between one representation and the other. Suppose you have a k*m-by-n matrix A stored in the first format above. We can rearrange it to move all common objects together (forming a stack of k m-by-n matrices), with:

A = [1 2 3; ...
     4 5 6; ...
     7 8 9; ...
    10 11 12; ...
    13 14 15; ...
    16 17 18; ...
    19 20 21; ...
    22 23 24];
n = 3;
B = reshape(permute(reshape(A,[],2,n),[3 2 1]),n,[])';

and the output is

B =

 1     2     3
13    14    15
 4     5     6
16    17    18
 7     8     9
19    20    21
10    11    12
22    23    24

To convert back:

C = reshape(permute(reshape(B',n,2,[]),[3 2 1]),[],n);

and we get again:

C =
 1     2     3
 4     5     6
 7     8     9
10    11    12
13    14    15
16    17    18
19    20    21
22    23    24

libigl Version 1.0 beta release

July 17th, 2014

libigl teaser

Daniele and I have completed a working tutorial for libigl. The course at SGP 2014 also encouraged us to advertise a proper Version 1.0 Beta release of libigl.

This change marks our confidence that libigl is ready for wide use in geometry processing research. While we still qualify our release with “Beta” to indicate that the library is a work in progress and documentation, for example, may not be as thorough as commercial or better established libraries, we’re continuing to develop libigl and encourage you to git pull our repo and stay up to date.

Here’s the table of contents to our tutorial. Each section comes with a working mini example main.cpp which demonstrates a specific function or routine available in libigl:

Small typo in “An Algorithm for the Construction of Intrinsic…”

July 11th, 2014

Upon implementing the Intrinsic Delaunay Laplacian according the helpful notes in “An Algorithm for the Construction of Intrinsic Delaunay Triangulations with Applications to Digital Geometry Processing” [Fisher et al 2006], I noticed a small typo. Specifically I’m following this version of the document. The issue is regarding the computation of the length of the “flipped edge” f. Currently the last step in section 2.2 uses the Law of Cosines, writing that:

f = sqrt(b^2 + c^2 - 2 cos(alpha + delta))

This should be:

f = sqrt(b^2 + c^2 - 2 b c cos(alpha + delta))

Pitfalls when implementing conformalized mean curvature flow

July 8th, 2014

Here’re some notes to avoid (common? obvious?) pitfalls when implementing conformalized mean curvature flow (“Can Mean-Curvature Flow Be Made Non-Singular?” [Kazhdan et al. 2012]).

Naively building a smoothing operation out of the identity matrix (data term) and a discrete pointwise Laplacian (smoothness term) would lead to iterations looking something like this:

L = cotmatrix(V,F);
I = speye(size(L));
U = V;
while true
  M = massmatrix(U,F,'barycentric');
  U = (I-delta*M\L)\U;
  tsurf(F,U)
  drawnow;
end

This is vary badly scaled as we need to first invert the mass matrix, introducing numerical problems if our mesh is irregular, and then solve the system.

Algebraically we can fix these issues by multiplying both sides of the system by our mass matrix:

L = cotmatrix(V,F);
U = V;
while true
  M = massmatrix(U,F,'barycentric');
  U = (M-delta*L)\(M*U);
  tsurf(F,U)
  drawnow;
end

This solution might converge for very small delta values, but eventually the surface positions U will converge on a single point and our floating point representation will run out of precision. In our context, the problem will be that entries in our diagonal mass matrix will degenerate because triangles in our mesh are degenerating.

So for numerical reasons we must normalize the positions at every step. There are various ways we could do this, here I show normalizing the surface to have unit surface area:

L = cotmatrix(V,F);
U = V;
while true
  M = massmatrix(U,F,'barycentric');
  U = (M-delta*L)\(M*U);
  U = U/sqrt(sum(doublearea(U,F))*0.5);
  tsurf(F,U)
  drawnow;
end

This is better, but still might only converge for small delta values and even then eventually diverge. The problem this time is a bit more subtle. Our mesh is unit area, but “floating” in space. Importantly, it’s floating away from the origin. Thus to store the coordinates, we’re wasting precision on each vertex by also storing this global translation. We eventually run into the same problem as before. Thus we should recenter the mesh around the origin at each step. Again, many ways to do this, but I show centering the centroid:

L = cotmatrix(V,F);
U = V;
while true
  M = massmatrix(U,F,'barycentric');
  U = (M-delta*L)\(M*U);
  U = U/sqrt(sum(doublearea(U,F))*0.5);
  U = bsxfun(@minus,U,centroid(U,F));
  tsurf(F,U)
  drawnow;
end

Depending on the quality of the original mesh we should now be able to take quite huge time steps and still see convergence (e.g. delta = 1e12!).

Important: Normalizing the surface area and recentering at the origin are necessary for numerics, not just visualization.

Note: Convergence means that our mesh flows to the sphere. However, the mesh may still slide around on the sphere. The problem being that any Möbius transformation is allowed by the flow. Eventually vertices may group up at a pole on the sphere and our mass matrix will degenerate.

Classic bump deformation example mesh

June 28th, 2014

This is code I end up re-writing every once in a while to demonstrate k-harmonic deformation on a simple 2d domain. The domain is just a square. This codes meshes the domain so that two inner circles show up in the edge set. Then I can select the vertices on and outside the outer circle and on and inside the inner circle to create a bump.

R = 1;
r = 0.15;
outer_sam = 200;
inner_sam = ceil(outer_sam/R*r);
inner_theta = linspace(0,2*pi,inner_sam+1)';
outer_theta = linspace(0,2*pi,outer_sam+1)';
P_inner = r*[cos(inner_theta(2:end)) sin(inner_theta(2:end))];
P_outer = R*[cos(outer_theta(2:end)) sin(outer_theta(2:end))];
E_inner = [1:size(P_inner,1);2:size(P_inner,1) 1]';
E_outer = [1:size(P_outer,1);2:size(P_outer,1) 1]';
P_bb = 1.1*[-R -R;-R R;R R;R -R];
E_bb = [1 2;2 3;3 4;4 1];
[P_bb,E_bb] = resample_polygon(P_bb,E_bb,2*pi*R/outer_sam);
P = [P_inner;P_outer;P_bb];
E = [E_inner;size(P_inner,1)+[E_outer;size(P_outer,1)+E_bb]];
[V,F] = triangle(P,E,[],'Flags',sprintf('-Yq32a%0.17f',(2*pi*R/outer_sam)^2),'Quiet'); 
tsurf(F,V)

bump domain matlab

Get polygon from user in matlab as list of points

June 25th, 2014

Simple as that:

P = getPosition(impoly);

Projecting out the null space of an underconstrained linear system

June 25th, 2014

This came up twice now recently, so I’m writing it out as a reference.

Suppose you have a system:

A x = b

where A is n by n, x is n by 1 and b is n by 1. However, rank(A) = m < n. This means there is not just one unique solution, but a family of solutions. Now, let’s assume we can build a basis for the null space of A. That is we can write an n by n-m matrix N such that

A N u = 0

for any n-m by 1 vector u. Then we can also write the matrix M = N which spans the space orthogonal to N. Then we can write any n by 1 vector in these bases:

x = N u + M v

Substituting this into our original system, we see that the choice of u doesn’t matter

A (N u + M v) = b
A M v = b

Now we have a system of n equations and m unknowns. This is over-constrained. Let’s choose the minimum in the least squares sense.

minimize |A M v - b|^2

Differentiating with respect to v and setting to zero we get

M' A' A M v = M' A' b

which we can solve for v.

Then we can recover x as

x = M v

We can show that indeed this satisfies the original equation by substituting this back in above:

A M v = b
A M ( M' A' A M ) \ M' A' b = b

Let Z = A M, then write this as

Z (Z' Z) \ Z b = b

But (Z’ Z) \ Z is simply the pseudo-inverse of Z so we know that Z (Z’ Z) \ Z = I and we have that b = b.

Eigen unaryExpr gotcha, need to set via assignment

June 23rd, 2014

I tried to use eigen’s unaryExpr with a function using references to alter a matrix:

MatrixXd X;
...
X.unaryExpr([](double & x){ x = 0});

But X was left unchanged. Thinking I might need to use a function with a return I tried:

X.unaryExpr([](double & x)->return{ return 0});

but still X was unchanged. It seems that unaryExpr is a const operation so I must assign the result back into X. This works:

X = X.unaryExpr([](double & x)->return{ return 0});

Continuously refresh markdown output

June 17th, 2014

I had a hard time finding a chrome plugin for a rendering markdown which supported LaTeX equations. Instead I’m using multimarkdown which supports equations with mathjax. However, to render the html output I need to open it in Chrome and refresh it when anything changes. So my pipeline looks like:

repeat
  edit .md files
  run multimarkdown command
  open/refresh .html file in chrome
end

I found this nice ruby/applescript script.
I’ve modified it to support markdown files and to run a command when a file changes:

#!/usr/bin/env ruby
# watch.rb by Brett Terpstra, 2011 <http://brettterpstra.com>
# with credit to Carlo Zottmann <https://github.com/carlo/haml-sass-file-watcher>

trap("SIGINT") { exit }

if ARGV.length < 2
  puts "Usage: #{$0} watch_folder keyword command"
  puts "Example: #{$0} . mywebproject"
  exit
end

dev_extension = 'dev'
filetypes = ['css','html','htm','php','rb','erb','less','js','md']
watch_folder = ARGV[0]
keyword = ARGV[1]
command = ARGV[2]
puts "Watching #{watch_folder} and subfolders for changes in project files..."

while true do
  first = true
  while true do
    files = []
    filetypes.each {|type|
      files += Dir.glob( File.join( watch_folder, "**", "*.#{type}" ) )
    }
    new_hash = files.collect {|f| [ f, File.stat(f).mtime.to_i ] }
    hash ||= new_hash
    diff_hash = new_hash - hash
    if not diff_hash.empty?
      if first and not command.empty?
        val=`#{command}`
        first = false
      else
        break
      end
    else
      sleep 1
    end
  end


  unless diff_hash.empty?
    hash = new_hash

    diff_hash.each do |df|
      puts "Detected change in #{df[0]}, refreshing"
      %x{osascript<<ENDGAME
            tell application "Google Chrome"
    set windowList to every window
    repeat with aWindow in windowList
        set tabList to every tab of aWindow
        repeat with atab in tabList
            if (URL of atab contains "#{keyword}") then
                tell atab to reload
            end if
        end repeat
    end repeat
        end tell
ENDGAME
}
    end
    sleep 1
  end

end

Save this in a file watch.rb and then call with something like:

watch.rb . readme "multimarkdown readme.md -o readme.html"