| brain activity log
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.
24.03.2006 Friday - diary
The GEOS system at Bassilichi is up, running and doing fine: another project
finished! It took several months of work in the end. The customer is satisfied as well as me.
Francesco already has another task waiting for me: the archive management system
must be rewritten from scratch (we're substituting an ancient unix-based system).
Neeeeeeeeeeext! :DD
Tried IE7 beta. It's a copy of firefox... well.. at least they have finally
discovered that tabbed browsing and transparent pngs are good (tm).
Comed back from work at 20, felt asleep for 5 minutes and woken up at 00 :D
Written a couple of blog lines and tried to fall asleep again. Failed.
Then decided to write down a template matrix class with a nice inversion
algorithm inside. Succeeded at nearly 6 am. Still failing to sleep watched
70% of "Taxi Driver", a cult movie. Finally felt asleep. Too good it's friday...
23.03.2006 Thursday - random thoughs about luck
The need for "fortunate randomness".
Well... sure: you have to help your luck, otherwise it will get to
nothing. But it's still randomness: the complex composition of events that lead
you to an advantageous situation. Often too complex to be described or even understood.
On the other hand... rationale.
In an example of good luck, a person winning a lottery would generally be considered lucky,
although a rationalist might point out that there was bound to be a winner sooner or later,
and there was actually nothing lucky about someone winning - it was merely a probabilistic event.
It is doubtful that the winner would agree with that analysis, however.
The belief in luck as a supernatural phenomenon is generally regarded by rationalists
as a form of magical thinking. However, there is evidence that people who believe
themselves to have good luck are more able to take advantage of fortunate chance events
in their lives, and to compensate for unfortunate chance events in their lives, than
people who believe that they have bad luck. This appears to be the result of positive
thinking altering their responses to these events.
If "good" and "bad" events occur at random to everyone, believers in good luck will
experience a net gain in their fortunes, and vice versa for believers in bad luck.
This is clearly likely to be self-reinforcing. Thus, although untrue, a belief in good
luck may actually be an adaptive meme.
"Luck is what happens when preparation meets opportunity" - Seneca, Roman Dramatist
"In my experience, there's no such thing as luck" - Obi-Wan Kenobi.
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.
| 