Inverse of a 3×3 block matrix

Recall the formula for the inverse of a 2×2 block matrix:

$$ \begin{aligned} \begin{bmatrix} A&B\\ C&D \end{bmatrix}^{-1} &= \begin{bmatrix} A^{-1}+A^{-1}BS^{-1}CA^{-1}&-A^{-1}BS^{-1}\\ -S^{-1}CA^{-1}&S^{-1} \end{bmatrix} \\ S &= D - C A^{-1} B \end{aligned} $$

Now consider a 3×3 block matrix

$$ \begin{aligned} X &= \begin{bmatrix} E&F&G\\ H&J&K\\ L&M&N \end{bmatrix} \end{aligned} $$

Apply the 2×2 block inverse formula, plugging in: $\tilde A=E$, $\tilde B=\begin{bmatrix}F&G\end{bmatrix}$, $\tilde C=\begin{bmatrix}H\\L\end{bmatrix}$, and $\tilde D=\begin{bmatrix}J&K\\M&N\end{bmatrix}$:

$$ \begin{aligned} X^{-1} &= \begin{bmatrix} E&\begin{bmatrix}F&G\end{bmatrix}\\ \begin{bmatrix}H\\L\end{bmatrix}&\begin{bmatrix}J&K\\M&N\end{bmatrix} \end{bmatrix}^{-1} \\&= \begin{bmatrix} E^{-1}+E^{-1}\begin{bmatrix}F&G\end{bmatrix}Z^{-1}\begin{bmatrix}H\\L\end{bmatrix}E^{-1}&-E^{-1}\begin{bmatrix}F&G\end{bmatrix}Z^{-1}\\ -Z^{-1}\begin{bmatrix}H\\L\end{bmatrix}E^{-1}&Z^{-1} \end{bmatrix} \\ Z &= \begin{bmatrix}J&K\\M&N\end{bmatrix} - \begin{bmatrix}H\\L\end{bmatrix} E^{-1} \begin{bmatrix}F&G\end{bmatrix} \end{aligned} $$

The factor $Z$, in turn, is another 2×2 block matrix.

$$ \begin{aligned} Z &= \begin{bmatrix}J&K\\M&N\end{bmatrix} - \begin{bmatrix}H\\L\end{bmatrix} E^{-1} \begin{bmatrix}F&G\end{bmatrix} \\ &= \begin{bmatrix}J&K\\M&N\end{bmatrix} - \begin{bmatrix} HE^{-1}F & HE^{-1}G \\ LE^{-1}F & LE^{-1}G \end{bmatrix} \\ &= \begin{bmatrix} J-HE^{-1}F & K-HE^{-1}G \\ M-LE^{-1}F & N-LE^{-1}G \end{bmatrix} \end{aligned} $$

Again apply again the 2×2 block inverse formula to get $Z^{-1}$, defining:

$$ \begin{aligned} A &= J-HE^{-1}F \\ B &= K-HE^{-1}G \\ C &= M-LE^{-1}F \\ D &= N-LE^{-1}G \end{aligned} $$$$ \begin{aligned} Z^{-1} &= \begin{bmatrix} A&B\\ C&D \end{bmatrix}^{-1} = \begin{bmatrix} A^{-1}+A^{-1}BS^{-1}CA^{-1}&-A^{-1}BS^{-1}\\ -S^{-1}CA^{-1}&S^{-1} \end{bmatrix} \\ S &= D - C A^{-1} B \end{aligned} $$

Further expanding the forumula for $X^{-1}$ in terms of this is tedious and somewhat unsatisfying. Define

$$ \begin{aligned} U &= G - FA^{-1}B \\ V &= L-CA^{-1}H \end{aligned} $$

then expand and simplify:

$$ \begin{aligned} E^{-1}&+E^{-1}\begin{bmatrix}F&G\end{bmatrix}Z^{-1}\begin{bmatrix}H\\L\end{bmatrix}E^{-1} \\&= E^{-1}\left\{I+ \begin{bmatrix}F&G\end{bmatrix} \begin{bmatrix}A^{-1}+A^{-1}BS^{-1}CA^{-1}&-A^{-1}BS^{-1}\\-S^{-1}CA^{-1}&S^{-1}\end{bmatrix} \begin{bmatrix}H\\L\end{bmatrix}E^{-1} \right\} \\&= E^{-1}\left\{I+ \begin{bmatrix}F&G\end{bmatrix} \begin{bmatrix}[A^{-1}+A^{-1}BS^{-1}CA^{-1}]H-A^{-1}BS^{-1}L\\-S^{-1}CA^{-1}H + S^{-1}L\end{bmatrix} E^{-1} \right\} \\&= E^{-1}\left\{I +\left\{ F\left[(A^{-1}+A^{-1}BS^{-1}CA^{-1})H-A^{-1}BS^{-1}L\right] +G\left[-S^{-1}CA^{-1}H + S^{-1}L\right] \right\}E^{-1} \right\} \\&= E^{-1}\left\{I + \left\{ FA^{-1}\left[H+BS^{-1}(CA^{-1}H-L)\right] -GS^{-1}[CA^{-1}H-L] \right\} E^{-1} \right\} \\&= E^{-1}\left\{I+\left\{FA^{-1}H+\left[FA^{-1}B-G\right]S^{-1}[CA^{-1}H-L]\right\}E^{-1}\right\} \\&= E^{-1}+E^{-1}\left\{FA^{-1}H+US^{-1}V\right\}E^{-1} \\ \\ -&E^{-1}\begin{bmatrix}F&G\end{bmatrix}Z^{-1} \\&= -E^{-1} \begin{bmatrix}F&G\end{bmatrix} \begin{bmatrix}A^{-1}+A^{-1}BS^{-1}CA^{-1}&-A^{-1}BS^{-1}\\-S^{-1}CA^{-1}&S^{-1}\end{bmatrix} \\&= -E^{-1}\begin{bmatrix} FA^{-1}+FA^{-1}BS^{-1}CA^{-1}-GS^{-1}CA^{-1} & -FA^{-1}BS^{-1}+GS^{-1} \end{bmatrix} \\&= \begin{bmatrix} -E^{-1}\left\{ F+(FA^{-1}B-G)S^{-1}C \right\} A^{-1} & E^{-1}[FA^{-1}B-G]S^{-1} \end{bmatrix} \\&= \begin{bmatrix} -E^{-1}\left[ F-US^{-1}C \right] A^{-1} & -E^{-1}US^{-1} \end{bmatrix} \\ \\ -&Z^{-1}\begin{bmatrix}H\\L\end{bmatrix}E^{-1} \\&=-\begin{bmatrix}A^{-1}+A^{-1}BS^{-1}CA^{-1}&-A^{-1}BS^{-1}\\-S^{-1}CA^{-1}&S^{-1}\end{bmatrix} \begin{bmatrix}H\\L\end{bmatrix}E^{-1} \\&= \begin{bmatrix} -A^{-1}H-A^{-1}BS^{-1}CA^{-1}H + A^{-1}BS^{-1}L \\ S^{-1}CA^{-1}H -S^{-1}L \end{bmatrix}E^{-1} \\&= \begin{bmatrix} -A^{-1}[H+BS^{-1}(CA^{-1}H-L)]E^{-1} \\ S^{-1}(CA^{-1}H-L)E^{-1} \end{bmatrix} \\&= \begin{bmatrix} -A^{-1}[H-BS^{-1}V]E^{-1} \\ -S^{-1}VE^{-1} \end{bmatrix} \end{aligned} $$

This gives the expanded formula:

$$ \begin{aligned} &X^{-1}= \\& \begin{bmatrix} E^{-1}+E^{-1}\left[FA^{-1}H+US^{-1}V\right]E^{-1} & -E^{-1}\left[F-US^{-1}C\right]A^{-1} & -E^{-1}US^{-1} \\ -A^{-1}[H-BS^{-1}V]E^{-1} & A^{-1}+A^{-1}BS^{-1}CA^{-1} & -A^{-1}BS^{-1} \\ -S^{-1}VE^{-1} & -S^{-1}CA^{-1} & S^{-1} \end{bmatrix} \end{aligned} $$


Random noise textures on the GPU in WebGL

Procedural rendering routines often need pseudorandom numbers. For graphics rendering, we usually prefer correct-ish and fast code over properly "correct" code. My approach to random number generation is usually to intermittently seed a weak pseudo-random number generator, like a linear-feedback shift register, from a stronger source like Javascript's Math.random(). 

On the GPU, things are not so simple: we want a pseudorandom algorithm that is fast, local, and parallel. For compatibility, we should also have this RNG use an 8-bit RGBA texture for its internal state. Bitops are a bit tricky in WebGL shaders (though not impossible), so we'll use a linear congruential generator.

The linear congruential generator is quite simple: it applies a linear transformation $ax+c$, then takes the remainder modulo some constant $m$. 

$$x_{n+1} = (a x_n + c) \operatorname{mod} m$$

Usually, we do this with integers, and $m$ is some prime number. While its possible to interpret the RGBA texture values as storing a byte [0,255] mapped to [0,1], running a weak RNG for each color channel is unsatisfactory. At best, the RNG has a period of 256, since the state stored between updates is effectively a single 8-bit unsigned integer for each color channel. 

Rather than try to interpret the RGBA vec4 as a 32-bit integer, I opted to mix state between adjacent pixels. This turns each channel in the texture into a giant RNG with as much state as you have pixels. I can't make any claims that this is optimal or correct, only that the results are fast and look nice. 

Initialize this texture with random values. Set the texture sampling to nearest-neighbor with periodic boundary conditions  (and no mipmaps, if it matters). In some instances, you may need to manually implement periodic boundary conditions for non power-of-two sized textures.

Here's the main code. It sums up pixels in the local neighborhood, multiplies them by some random float, then stores them back in the texture modulo 1.0 (that's modulo 256 if we're thinking of the color data as unsigned integers). I forget how I chose the constant for the linear recurrence, but it was almost surely a variant of the "try random things, smash the keyboard, and see what works" algorithm.

void main() {
    vec2 scale = 1./vec2(W,H);
    vec4 c0 = texture2D(data,(XY+vec2( 0, 0))*scale);
    vec4 c1 = texture2D(data,(XY+vec2( 0, 1))*scale);
    vec4 c2 = texture2D(data,(XY+vec2( 1, 0))*scale);
    vec4 c3 = texture2D(data,(XY+vec2(-1, 0))*scale);
    vec4 c4 = texture2D(data,(XY+vec2( 0,-1))*scale);
    vec3 color = mod((c0+c1+c2+c3+c4).rgb
    gl_FragColor = vec4(color,1.);

This gives you a uniform RNG for each color channel in the texture to play with. If you want Gaussian distributed data, you can use the Box-Muller transform to convert the (R,G) and (B,A) channels in a pair of Gaussian random numbers each.  

This random noise shader is up as an example on the webgpgpu repository. Click on the image below to view an example of Gaussian color noise: 


WebGL shader rules for rendering dendriform growth

I've been playing around with cellular automata for procedural generation in WebGL. Here's an algorithm for rendering dendriform growth. I was aiming for rivers, but it came out like neurons.

[view in browser]

It's built atop a diffusion-limited aggregation cellular automata. Algorithm sketch: 

  • Each cell looks at the 4 nearest neighbors on the grid. We also check the 4 corners ("far" neighbors) to avoid creating loops.
  • Random numbers are provided by a separate noise texture
  • State is stored in the channels of an RGBA texture. Each [0,1] channel maps to a [0,255] byte value
  • The texture is seeded with at least one "active" pixel (existing water/lake/cell body). The seed is set to active (green=1.0)
  • Green channel runs a a diffusion limited aggregation: 
    • If exactly one nearby pixel is active, then we have a 5% of becoming active. 
    • Avoid creating loops: don't activate a pixel if more than one "far" neighbor is already active.
    • Because pixels act in parallel and at random, sometimes two pixels turn on in the same iteration, creating a loop. Detect and remove these after-the-fact, if they happen. 
  • Keep track of distance from the "seed" region in blue channel
    • Seeds are initialized with blue channel = 255
    • Newly activated cells get blue-1
    • This counts down to zero, at which point expansion stops
  • The texture is also initialized with various water sources / "synapses" scattered about it
    • This is stored in the alpha channel
    • If an active cell hits a source, it becomes a "river"
    • If a cell is active, and an adjacent cell is a "river", and is also further away from the source (check distance in blue channel), then this cell also becomes a "river"
  • A separate shader checks which cells are rivers, and generates tile ID values which are then passed to a tile shader.

Pixel-art-style Conway's game-of-life using a tile shader in WebGL

Implementation of Conway's game of life cellular automata with a retro-style shader. Added as example on webgpgpu Github.

Sketch of WebGL implementation:

  • Underlying cellular automata runs Conway's game of life
  • States as game tiles: new cells: small blue; live cells: blue; died: skulls
  • Life diffuses out coloring nearby areas
  • A single RGBA texture stores the game state. We render to/from this texture as a framebuffer.
  • Each color channel in [0,1] can be mapped to a [0,255] byte value
  • The green channel stores the game of life state
  • The blue channel stores the diffusing field
  • The red channel outputs a tile ID value that is passed to a tile shader to render the game

Retro-styled forest-fire cellular automata in WebGL

This is a cellular automata analogue of the model that we used for developmental retinal waves in this paper

Key rules

  • Vegetation "grows" occasionally (flip a coin to increment amount of trees)
  • Fires have a small chance of starting spontaneously
  • Intensity of ignited fires decay down to zero, at which point they leave ash
  • Fires spread to adjacent cells with a probability that depends on their intensity
  • Probability of a cell igniting given a fire is related to how much vegetation there is

Sketch of WebGL implementation:

  • A single RGBA texture stores the game state
  • Each color channel in [0,1] can be mapped to a [0,255] byte value
  • A "noise" texture creates pseudorandom numbers
  • The green channel stores the level of vegetation (BURNT,NONE,GRASS,SCRUB,TREE)
  • The blue channel stores the intensity of the fire
  • The red channel outputs a tile ID value that is passed to a tile shader to render the game


Algorithms for rendering quasicrystal tilings

Quasicrystals are aperiodic spatial tilings, and 2D quasicrystals have been used for centuries in Islamic art. I find these patterns beautiful, and I think more people should know how to render them.

You can render 2D quasicrystals as a sum of plane waves. They can also be constructed using aperiodic tilings (e.g. 4-fold Ammann-Beenker tiles and 5-fold Penrose Tiles), and there are fractal constructions that recursively subdivide a set of tiles using a smaller set of the same tiles (e.g. 5- and 7-fold tilings). Many more constructions can be found on the Tilings Encyclopedia

This post focuses on the "cut and project" construction, which is easy to implement in code. The routines shown here can be found in this iPython notebook on Github.

The cut-and-project construction

The simplest way to render crisp, geometric quasicrystals is the "cut and project" construction. This views the quasicrystal as a two-dimensional projection of a planar slice through a N-dimensional hyper-lattice. Let's implement this algorithm. 

A (hyper) lattice

First, we need to define the idea of a hyperlattice. This is a N-dimensional generalization of a cubic crystal lattice. It tiles ND space using ND hypercubes. We can represent the center of each hypercube (hypercell? hyperpixel?) as a N-D vector of integers
\mathbf q_0 = [ x_1, ... x_N ],\,\,x_i\in\mathbb Z
The associated N-cube, then, includes all points on a unit hypercube $\left[-\tfrac 1 2, \tfrac 1 2\right]^N$ offset to location $\mathbf q_0$:
\mathbf q = \left[-\tfrac 1 2, \tfrac 1 2\right]^N + \mathbf q_0

An irrational projection onto a 2D plane

To construct the quasicrystal, we cut through this ND lattice with a 2D plane. To get a aperiodic slice, this plane must not align with any periodic direction in the hyperlattice. 

A good choice is to project each direction of the hyperlattice onto an evenly-spaced set of basis directions, centered at the origin. We define a 2D $(s,t)$ plane $\mathcal P$ in terms of two N-D basis vectors, $\mathbf u, \mathbf v \in \mathbb R^N$.

\mathbf p &= \mathbf u \cdot s + \mathbf v \cdot t,\,\,s,t\in\mathbb R
\mathbf u &= [ \cos(\theta_0), \dots , \cos(\theta_{N-1}) ]
\mathbf v &= [ \sin(\theta_0), \dots , \sin(\theta_{N-1}) ]
\theta_i &= \frac {\pi\cdot i}{N}


Inverting the "Simple Sabotage Field Manual"

In 1944, the US Office of Strategic Services published an internal document: the "Simple Sabotage Field Manual". 

In addition to specific techniques for breaking machines, it included suggestions for how to cripple an organization. 

All of the suggestions are things that already happen spontaneously, because organizing large numbers of people reliably and safely is hard.

Rather than re-state the "sabotage" tips, I thought I'd invert each one to create guidelines for positive change in an organization:

  • Always share your competencies and knowledge freely and liberally.
  • Expedite when it is safe to do so, and whenever failing to do so would lead to something worse.
  • Rules that don't serve a purpose can be changed, or broken in emergencies (but this must be done with extreme thought, caution, and community support, see Chesterton's fence).
  • No speeches. Be brief and clear, and say nothing without purpose.
  • No committees, but if you must, never more than four persons.
  • Stay on mission. Do not re-open closed (or no longer relevant) matters.
  • With regards to language: be liberal in what you accept and conservative in what you give out; if the words communicated understanding, spend no time revising them.
  • Leadership requires action; Take thoughtful and calculated risks (hedging if appropriate). Beware inertia and cowardice in the name of caution.
  • Prioritize important jobs first, and assign them to the most competent and efficient.
  • Avoid perfectionism, be pragmatic.
  • Allocate promotions and responsibilities according to ability (but be kind and take care of all).
  • Conferences are a waste of time, if you already know what you need to do.
  • Minimize admin bottlenecks, but If a process really needs N persons to sign off, build in redundancy: ensure that at least 2N persons are trained and authorized to do so.
  • Be diligent and serious in your work.
  • Avoid distractions, minimize interruptions from co-wokers and avoid breaking up the workday.
  • Learn how to build and maintain good tools, and how to do good work with bad tools.