Robin

Matrices and Determinants

Posted in Robin on September 2nd, 2010 by Robin – Be the first to comment

Let denote an invertible matrix, and let denote its inverse.
Suppose that we are given a column vector and that we wish to solve the equation for the unknown vector .

Since is invertible with inverse , this is equivalent to:


In co-ordinates:

This would be great if we had a nice expression for the inverse of a matrix.

Now, let denote the columns of .
That is:


If denotes the standard basis then:

That is, the columns of the matrix are the images of the standard basis vectors.

Another way of looking at the problem is as follows. By definition, we have:


To solve the equation for we must find unknowns such that:

In other words, we are trying to find the co-ordinates of the vector in terms of a new basis.

Recalling now that the determinant is multi-linear and anti-symmetric, consider the expression:



Let us now rewrite the equation as the following “formal” determinant , which has entries in the first row the standard basis vectors, and additional “projective symbol” :

Expanding along the top row, whilst making use of the previous remark (and being careful with the signs), we get:




which after “deprojectivizing” gives us exactly the vector which we were searching for.

Next, by expanding along the -th column, we obtain:


where denotes the determinant of the minor obtained from by removing the -th row and the -th column. From which we recover the well-known formula for the inverse of a matrix:

Special Relativity and Hyperbolae

Posted in Robin on August 27th, 2010 by Robin – Be the first to comment

The matrix for the Lorenz Transformation, in one space dimension and one time dimension (with ) is given by:

I shan’t derive this matrix from the “principle of relativity” which states that the laws of physics are identical in all inertial frames of reference, for this has been done very well elsewhere, and is not the purpose of this post.

Let us return instead to our good friend, the two by two rotation matrix, in order to make a few observations:

Consider a right angled triangle with sides lengths , and . If you remember the definition of your basic trigonometric functions, as well as Pythagorus’ theorem, then you will “see” that we have:


This gives us the following alternative parameterization for the rotation matrices, where :


It looks rather like the Lorenz transformation, doesn’t it? Though not quite.

If you recall now the definition of the hyperbolic trigonometric functions, you will “see” that and , and as a consequence, we have:


This in turn gives us an alternative parameterization for the Lorenz Transformation, where :

The upshot of this is that if you are willing to work over rather than over , and to think of your space dimension as purely real, and your time dimension as purely imaginary, then applying a “Lorenz boost” to your frame of reference is equivalent to rotating through an imaginary angle.


In the same way that a regular rotation maps circles to circles in two space dimensions, a “hyperbolic rotation” sends hyperbolae to hyperbolae in one space dimension and one time dimension. The following mathematica code will let you play around with this idea.


points = {{1, 0}, {2, 0}, {3, 0}, {4, 0}};

A[theta_] := {{Cos[theta], -Sin[theta]}, {Sin[theta], Cos[theta]}};
B[v_] := {{1/Sqrt[1 - v^2], v/Sqrt[1 - v^2]}, {v/Sqrt[1 - v^2],
    1/Sqrt[1 - v^2]}};

Plot[y /. Table[Solve[x^2 + y^2 == k^2, y], {k, 5}], {x, -5, 5},
  PlotRange -> {-5, 5}, AspectRatio -> 1];
Manipulate[
 Show[%, ListPlot[points . A[theta], PlotStyle -> PointSize[0.025],
   PlotRange -> {{-5, 5}, {-5, 5}}, AxesOrigin -> {0, 0},
   AspectRatio -> 1]], {{theta, 0}, -Pi, Pi}]

Plot[y /. Table[Solve[x^2 - y^2 == k^2, y], {k, 0, 8}], {x, -10, 10},
  PlotRange -> {-10, 10}, AspectRatio -> 1];
Manipulate[
 Show[%, ListPlot[points . B[theta], PlotStyle -> PointSize[0.025],
   PlotRange -> {{-5, 5}, {-5, 5}}, AxesOrigin -> {0, 0},
   AspectRatio -> 1]], {{theta, 0}, -0.9, 0.9}]

Since a particle moving at the speed of light is represented in this picture by the degenerate hyperbola which form the asymptotes of the rest of the family, it is clear now why the speed of light is constant in all inertial frames of reference.

But what’s really going on here? Let be the matrix representing a Minkowski bilinear form with one space dimension and one time dimension:


Then any matrix of the form:

will have the property that:

Compare this to the Euclidean case, in which the bilinear form is represented by the identity matrix. Now we have that any matrix of the form:


has the property that:

In co-ordinates, if and represent two different particles in two dimensional space, then the Euclidean inner product is usually represented as:

In the Minkowski picture, where time is imaginary, a “point” represents a particle located at a given position in one dimensional space, at a given time. The “inner product” of two “points” and is given by:

Conic Sections and Projective Space

Posted in Robin on August 20th, 2010 by Robin – Be the first to comment

I’m sure you’ve all seen the equation of a circle before:


You’ve probably also seen the equation of a hyperbola:

Or perhaps you’re more familliar with seeing the asymptotes at 45 degree angles to the co-ordinate axes, as in:

All these are special cases of the most general quadratic equation, which can be written in the following projective form:

Note that the by matrix in the middle is symmetric. As an example, the Family of circles centered at the origin, with radius are associated to the family of matrices of the form:

The first form of the hyperbola, with asymptotes parallel to the co-ordinate axes, corresponds to the matrix:


while the second hyperbola, with asymptotes at 45 degrees, corresponds to the matrix:

More generally, the matrix associated to a hyperbola whose axes have been tilded by an angle of clockwise from the diagonal is:



The following mathematica code will let you play around with these rotations.

   v = Transpose[{{x, y, 1}}];
   similtude[A_, M_] := Transpose[M].A.M
   equation[A_] := Solve[similtude[A, v][[1]][[1]] == 0, y]

   R[theta_] := {{Cos[theta], -Sin[theta], 0},
                 {Sin[theta], Cos[theta], 0},
                 {0, 0, 1}};
   H = {{1,0,0},{0,-1,0},{0,0,-1}}
   W[theta_] := similtude[H, R[theta]]

   Manipulate[
     Plot[Evaluate[y /. equation[W[theta]]], {x, -4, 4},
       PlotRange -> {{-4, 4}, {-4, 4}},
       AspectRatio -> 1],
   {{theta, 0}, -Pi, Pi}]

A few remarks about the code. In Mathematica square brackets are reserved for passing arguments into functions, curly brackets for defining lists, and normal round brackets for grouping terms in algebraic expressions.

Double square backets are used for indexing elements of a list. Matrices are represented as lists of list. The first list is the first row, the second list is the second row, etc…

The matrix multiplication operator is simply a dot. When defining your own functions, the variables which can be passed into a function must be followed by an underscore.

The built in functions, Transpose and Solve do more or less what you would expect. The only subtlety is that the output of Solve is a transformation rule. The /. notation replaces the left hand side with the right hand side in the transformation rule.

The built in function Plot is responnsible for drawing static plots. Leaving a parameter free (in this case the angle theta) and wrapping it inside the Manipulate command allows you to vary the parameter with the mouse and view the resulting animation.

Recall that translation can be thought of an operation in projective space represented by the matrix:

Now, starting from the matrix for a circle of radius one, centered at the origin:


and applying the following similarity transform:



We obtain the matrix for a circle of radius one, centered at the point . The following code will let you play around with this:

  T[a_, b_] = {{1, 0, -a}, {0, 1, -b}, {0, 0, 1}};
  CC = {{1,0,0},{0,1,0},{0,0,-1}}
  M[a_, b_] := similtude[C, T[a, b]]
  Manipulate[
    Plot[Evaluate[y /. equation[CC]], {x, -2, 2},
      PlotRange -> {{-2, 2}, {-2, 2}}, AspectRatio -> 1],
  {{a, 0}, -1, 1}, {{b, 0}, -1, 1}]

Now for the fun part. Let us define the following “projective rotation”

Using this, we can smoothly deform a circle into a hyperbola. Check it out:

   P[theta_] := {{Cos[theta], 0, -Sin[theta]}, {0, 1, 0},
                 {Sin[theta], 0, Cos[theta]}};
   S[theta_] := similtude[CC, P[theta]]

   Manipulate[
     Plot[Evaluate[y /. equation[S[theta]]], {x, -10, 10},
       PlotRange -> {{-10, 10}, {-10, 10}},
       AspectRatio -> 1],
   {{theta, 0}, -Pi, Pi}]

But what’s really going on here geometrically? The following three dimensional quadratic equation describes a mathematical cone (which is more like a pair of icecream cones, touching at their tips):

By setting we are effectively taking the intersection with a plane. As you may know, all quadrics can be obtained by intersecting a cone with some given plane (hence the term conic section). The plane gives a circle, while the plane gives a hyperbola. Our “projective rotation” smoothly rotates the axis to the axis.

Translations as Projective Transformations

Posted in Robin on August 11th, 2010 by Robin – Be the first to comment

You probably already know that in two dimensional space, the geometric act of rotating a vector through a given angle can be represented algebraically by the matrix:

You might have also asked yourself “What about translation?” That’s a nice geometric operation, how can we represent it algebraically in terms of matrices?

The problem with translation is that it doesn’t preserve the origin, and is thus not a linear transformation. It is however a projective transformation, and so by working in projective space it is still possible to think of translations in terms of matrix operations. This is a common trick used in computer graphics packages such as openGL.

There are many different ways to think about two dimensional projective space. Geometrically it is the set of lines in three dimensional space which pass through the origin. Topologically its the two dimensional sphere with antipodal points identified. Think of the fact that a line through the origin in intersects the unit sphere in exactly two points.

Algebraically, it is the set of triples under the equivalence relation for all .

Another useful way to think about projective two-space is as regular two dimensional space with an additional “circle* at infinity”, that is an additional point at infinity for every possible direction which you might head out in from the origin. Imagine the hyperplane sitting inside and notice how any line passing through the origin, which does not lie in the two dimensional subspace spanned by the and directions, intersects this plane in exactly one point.

Working with the algebraic definition, a two dimensional projective transformation is just a by matrix acting on projective co-ordinates. The equivalence relation on the underlying projective space implies that two matrices which differ only by a scalar multiple represent the same projective transformation.

Suppose that we want to send the arbitrary point to the new point which is shifted across two steps and down three steps.

The first step is to embed the two dimensional plane, which we are really interested in, into the somewhat more abstract and difficult to understand two dimensional projective space, where we are going to perform our transformation. This is accomplished by the map:

The “inverse” of this embedding is the map:


I say “inverse” in inverted commas, because this map is not defined everywhere. In particular its not defined when .

Alternatively, if you have trouble visualizing projective space, just imagine embedding into as the hyperplane given by the equation . The only thing we’re really missing in this picture is the “circle at infinity” which will be important later when we consider quadrics, but does not really matter just now while we are only concerned with translation operators (which send infinity to infinity anyhow).

Consider now the following projective transformation:


(or if you prefer, just think of it as a regular three dimensional transformation, which happens to map the hyperplane back to itself)

If we apply this to the image of our point under our funny embedding, then we get:

Et Voila! After taking the inverse of the embedding, this is exactly the image of the point that we were looking for.

* As correctly pointed out by Cale Gibbard, its not actually a circle at infinity, but a projective line. If two people head out in opposite directions along the plane they will reach the same point at infinity.

Back from the Dead

Posted in Robin on August 5th, 2010 by Robin – Be the first to comment

This blog has been dead for some time, but since I am now officially on the job market, I have decided to re-open it. I plan to post something new at least once a week.

Data Mining

Posted in Links, Robin on January 28th, 2010 by Robin – Be the first to comment

For anybody who is interested, here are two online video lecture series from Stanford University:

Statistical Aspects of Data Mining

Machine Learning

Also, the data set for the netflix prize.

Why’s (poignant) guide to Ruby

Posted in Links, Robin on January 21st, 2010 by Robin – Be the first to comment

Despite numerous other unfinished projects, I have decided to learn Ruby.

Regardless of whether you have any actual interest in learning Ruby or not, if your laughing muscles need a work out, may I recommend:

Why’s (poignant) guide to Ruby

It even has a soundtrack.

Another nice analogy

Posted in Links, Robin on January 7th, 2010 by Robin – Be the first to comment

The field with one element

A nice analogy

Posted in Links, Robin on January 6th, 2010 by Robin – Be the first to comment

The Anatomy of Integers and Permutations.

The best part is – they even made it into a screen play.

Matrices in Haskell II

Posted in Robin on December 14th, 2009 by Robin – Be the first to comment

What’s the use of a matrix if you can’t find its characteristic polynomial? So, a very primitive haskell type for polynomials:

 
import MyPermutations
import MyMatrices
 
type Polynomial = [Int]
 
polySum :: Polynomial -> Polynomial -> Polynomial
polySum p q = if (length p >= length q) 
	then zipWith (+) p (q ++ [0,0..])
	else zipWith (+) q (p ++ [0,0..])
 
polyProduct :: Polynomial -> Polynomial -> Polynomial
polyProduct p q = map (polyProduct' p q) [0..((length p) + (length q) - 2)]
 
polyProduct' :: Polynomial -> Polynomial -> Int -> Int
polyProduct' p q k = sum [ ((p++[0,0..])!!i) * ((q++[0,0..])!!(k-i)) | 
                                         i <- [0..k] ]
 
.

Its not too much effort to generalize matrices with integer entries to matrices with polynomial entries:

 
type MatrixPoly = [[Polynomial]]
 
polyMatTranspose :: MatrixPoly -> MatrixPoly
polyMatTranspose ([]:xs) = []
polyMatTranspose xs = (map head xs):(polyMatTranspose (map tail xs))
 
polyMatProduct :: MatrixPoly -> MatrixPoly -> MatrixPoly
polyMatProduct m1 m2 = polyMatProduct' m1 (polyMatTranspose m2)
 
polyMatProduct' :: MatrixPoly -> MatrixPoly -> MatrixPoly
polyMatProduct' [] m = []
polyMatProduct' (r:rs) m = (polyFindRow' r m):(polyMatProduct' rs m)
 
polyFindRow' :: [Polynomial] -> MatrixPoly -> [Polynomial]
polyFindRow' r m = map (foldr (polySum) [0]) (map (zipWith (polyProduct) r) m)
 
polyDeterminant :: MatrixPoly -> Polynomial
polyDeterminant m = foldr (polySum) [0] (map (polyDetTerm m) 
                                           (allPermutedLists (length m)))
 
polyDetTerm :: MatrixPoly -> PermutedList -> Polynomial
polyDetTerm m p = map (\x -> (-1)^(sign p) * x) (polyDetTerm' m p)
 
polyDetTerm' :: MatrixPoly -> PermutedList -> Polynomial
polyDetTerm' _ [] = [1]
polyDetTerm' (m:ms) (p:ps) = polyProduct (m!!(p-1)) (polyDetTerm' ms ps)
 
polyMinor :: MatrixPoly -> [Int] -> [Int] -> MatrixPoly
polyMinor m rows cols = [ [ (m!!(i-1))!!(j-1) | j <- cols ] | i <- rows ] 
 
almostPolyMatInverse :: MatrixPoly -> MatrixPoly
almostPolyMatInverse m =  [ [ polyDeterminant $ polyMinor m (foo n i) (foo n j) 
                                        |  i <- [1..n] ] | j <- [1..n]]
	where n = length m
 
.

Finally:

 
charpoly :: Matrix -> Polynomial
charpoly m = polyDeterminant $ matPolyAdd (matToPolyMat m) 
      (diagonalPolyMat (take (length m) (repeat [0,-1])))
 
matToPolyMat :: Matrix -> MatrixPoly
matToPolyMat = map (map (\x -> [x])) 
 
matPolyAdd :: MatrixPoly -> MatrixPoly -> MatrixPoly
matPolyAdd = zipWith (zipWith (polySum))
 
diagonalPolyMat :: [Polynomial] -> MatrixPoly
diagonalPolyMat xs = diagonalPolyMat' (length xs) 1 xs
 
diagonalPolyMat' :: Int -> Int -> [Polynomial] -> MatrixPoly
diagonalPolyMat' _ _ [] = []
diagonalPolyMat' n k (x:xs) = ((take (k-1) (repeat [0])) ++ [x] 
         ++ (take (n-k) (repeat [0]))):(diagonalPolyMat' n (k+1) xs)
 
.