/ home / blog
brain activity log

01.04.2006 Saturday - April Fools

1st April. In many countries this is the day of "Jokes". "Prima Aprilis: uwazaj bo sie pomylisz".. which is in Polish and means "1 April, be careful because you might be wrong". "April Fools" and "Pesce d'Aprile". The (net-)jokes already started.

Heh.. after you read some of them, all the of the web starts to look like a joke :D

An Insightful post on slashdot stated:

Try date -u and you will see that it is in fact April Fool's Day.

Note that April Fool's Day, as defined by the International April Pranksters Association, goes by UTC, not by local time zones. Because IAPA is not widely recognized as an international standards body and many people's problems to understand time zones, this has led to some problems since its introduction in 2002. It is especially uncertain for regional publications, which are reluctant to adopt the new standard because they fear to irritate their audience. One example is the Hubsborough Gazette, which famously spread confusion on the evening of March 31st, 2004 (EST) when an article claiming that aliens have attacked the Whitehouse appeared on their website. Despite the seemingly obvious nature of the hoax, many believed it and called the authorities or local clergery for guidance. One family even is reported to have spent two weeks in their backyard bunker. Since then, the editor has announced that they will only publish April Fool's articles during the hours when April 1st of their local time zone and UTC overlap, and take down articles afterwards. Many publications have followed their example in the following years.

The guys from Centrica fooled me about bex2 emitting all zero size output... that would be a problem. It is a 1st April joke..isn't it ? :D

31.03.2006 Friday - diary and other goodies

Hum.. a lot of images (in random order):

  • The Diophantine Equation
  • A nearly black cat
  • The Youla-Kucera theorem (kinda hard proof for the multivariable case)
  • A pair of sweet lips :)
  • The missed Frantz sister's party: mea culpa :(
  • Francesca in the lecture room
  • Silvia pressing for messaging
  • Samuel spending up to the last cent for a drum set
  • The two girls near Skile tripping three times around the block
  • A sensation close to deja-vu
  • Iakko laughing hard after drinking 4 Hoegarden glasses :D
  • ...
Nice four days :)

The promised LU decomposition code is here. Pretty straightforward once the formula is known. The file contains also the algorithms for matrix determinant (computed either by LU or by a recursive adaptation of the Liebnitz fomula) and for matrix inversion. Enjoy :)

"I have a monster in me and if I wake him up he is probably going to kill me"

27.03.2006 Monday - The LU decomposition of a matrix

Back to my fast determinant computed by using the LU decomposition of a matrix.

The idea is to decompose the matrix A in a product of two matrixes L and U in that L is lower triangular and U is upper triangular. We then know that det(A) = det(L) * det(U) and the determinants of triangular matrixes are trivial to compute. The LU decomposition is not unique so we can finally choose some arbitrary elements of L or U. We choose the diagonal of L to be made of all ones so the determinant of L is directly computed as 1.

Determining L and U such that A = LU is a system of m^2 equations (where m is the rank of A) in 2 * m^2 unknown variables. The bounds of L being lower triangular with a diagonal of all ones and U being upper triangular reduce the number of unknown variables to m^2. The system looks like the following one (3 x 3 case):

Which can be expanded to:

This shows us that we can compute each row of L and U independently by solving a m-equantion in m-variables system. In the 3x3 case the three systems look like:

Which can be solved by inverting the coefficient matrix or by Gauss-Jordan manipulation. The Gauss-Jordan solution for the 3x3 case looks like:

There's obviously a pattern in the solution vector and it can be written down as follows:

Now we have a LU decomposition and can compute the determinant just by multiplying together the diagonal elements of U. Nice eh ?

This is pretty fast. While the Liebnitz method costs like n!, the LU method costs like n^3 which is faaaaaar bettter :)

Now one should prove that the LU decomposition exists for every non-singular matrix and verify the numeric stability of the method... I'll do it in the next days while implementing this stuff and report you later.

26.03.2006 Sunday - 1 is a number

  • 1 is a number
  • Every number has a successor (next one) that itself is a number

  • The number "before" the successor is the predecessor (of the successor)
This allows us to define the set N of natural numbers. We give each number a "name" (one, two, three...) that is defined by rules beyond the scope of this document.

  • the predecessor of 1 is 0
We include 0 in the set N to allow the following definitions.
  • We define the operation of sum as follows:
    The sum of number a and number b is computed by applying the "successor" rule to a and the "predecessor" rule to b until b reaches 0.
For example summing four to three leads to the following steps:
  • successor of four is five, predecessor of three is two
  • successor of five is six, predecessor of two is one
  • successor of six is seven, predecessor of one is zero
  • the result is seven

  • We define the operation of subtraction as the inverse of sum. Being pedantic one can define it in terms of successor and predecessor itself: subtracting b from a means applying the predecessor rule to both a and b until b reaches zero.
Oops... we can't subtract b from a when b is in the successor chain of a: we haven't defined the predecessor of zero. We define it now as well as all the other predecessors. This leads us to the definition of the set Z of relative integers. Z obviously includes N.

  • We define multiplication of two numbers a and b as summing b to zero a times. Being pedantic: we start with an accumulator of 0 and sum b to it while applying the predecessor rule to a until it reaches zero.
  • We define division as the inverse operation of multiplication. If a multiplied b given c then c divided b gives a and c divided a gives b.
Oops... we can now attempt to divide b by a when there is no number in Z the multiplied a gives b.
  • Such a "thing" is still a number and is exactly b/a: a rational number.
This leads us to the definition of the set Q of rational numbers: the ones that can be rappresented by a fraction (ratio). Q obviously includes Z and thus N.

  • We define the "power" operation of numbers a and b as muliplying 1 for a b times. Thus a to the power of b is computed as multiplying 1 for a and applying the predecessor rule to b until zero is reached.
  • We define the extracion of the N-th root as the inverse of the power operation. If a to the power of b given c then the b-th root of c gives a.
Oops... the N-th root of certain numbers is not inside Q: you can attempt to extract the N-th root of a number x when there is no number y in Q that leads y^n=x. The square root of 2 is such a number. This leads us to the definition of the set R of real numbers. Obviously R contains Q and thus Z and N. R is dense: between two numbers in R there is always another number.

Oops... this still does not allow the extraction of certain N-th roots of negative numbers. The square root of -1 doesn't lie in R. In other words: there is no such number x in R that x^2=-1. We define such a number as the immaginary unity j (engineers use j :) and extend the set R with it. This leads us to the definition of the set C of complex numbers: the ones that have a real part (that lies in R) and an immaginary part (that doesn't lie in R). We rappresent such numbers as a+jb.


Just to make sure you know :)

Now go and conquer the world!

25.03.2006 Saturday - matrixes

Yesterday and today my fundamental problem was inverting a matrix and computing its determinant. I have written a nice template matrix class (I'll publish it on this site sooner or later).

Computing the inverse of a matrix is pretty straightforward. It can be done by Gauss elimination and has computational cost comparable to n^3 where n is the rank. Writing the Gauss elimination algo took a couple of hours yesterday night but afterwards I was able to invert a 1000x1000 matrix in matters of seconds (Athlon64 3500).

A huger problem is the one of computing the determinant. The Liebnitz formula is too complex to implement since it involves enumerating all the permutations of the matrix columns. Enumerations of n elements are n! and such complexity is obviously out of range.

Using the recursive determinant formula is also a heavy approach since it requires a lot of time and a really huge amounts of memory. You have to compute ALL the minors of the matrix... which is too much.

I've written down a hybrid method that uses a sort of Liebnitz approach with the first row and a recursive method for the minors. It's still too much. Computing a 20x20 determinant is not feasible in terms of time (it requires several hours) and it also eats a lot of memory since you need to hold 20x20 + 19x19 + 18x18 + 17x17 + ... + 1x1 matrixes in memory at once.

I've readed of the LU decomposition approach which should be comparable to the inversion in terms of complexity and memory usage. The idea is to decompose our matrix in a product of two matrices: L and U. L should be lower triangular and U upper triangular. Once such a decomposition is known the determinant of the matrix is det(L) * del(U) which are both products of the elements on the diagonal. There is an approach in that L has all 1 elements on the diagonal so det(L) is even 1 (no need to compute it at all).

The problem, now, is to find the LU decomposition. I still have to study better the algorithm. This will be a task for tomorrow tough. Tonight I'll also ask Valeria about it. She's a matematician so she probably will be able to give me some hints.

want more ?

... really ? :D

Browse around then.

You're viewing 5 posts per page: you can view more or less, if you want.

The entries marked in red are the ones you're viewing now.