Matrix equation

Jani wrote:
>> ...
>> 1. R_k and R_k^-1 are circulant kxk matrices (because R_k is one, and
>> R_k is positive definite => invertible)
>
>Unfortunately they are not. Altough I'm willing to do that approximation
>(because R is "almost circulant"), if it considerably helps in this
>situation.

Hmm, that makes it more difficult (perhaps). I'm not sure how one can
generally find eigendecompositions of Toeplitz matrices (R is normal,
so there exists a vector space basis of eigenvectors of R). Once you
have the general eigendecomposition of R_k you can perhaps also find it
for R_k^-1 - S_k (at least finding it for R_k^-1 is easy). If this
operator has only one non-zero eigenvalue, all you have to do is
project r onto that eigenspace and multiply the projection with the
corresponding eigenvalue to get e.

For example, if R is a diagonal matrix, ie. R_k = c 1_k, then R_k^-1 -
S_k has eigenvalue 1/c with eigenspace spanned by e_k, the k-th
standard normal basis vector.

Perhaps you can find a proof by induction: show that it (the fact that
R_k^-1 - S_k has only one non-zero eigenvalue)  holds for k=2, where it
shouldn't be too difficult to compute everything explicitely, and then
proceed with k->k+1.

FWIW,
Andor

Jani Huhtanen wrote:
> abariska@student.ethz.ch wrote:
> 
>> ...
>> 1. R_k and R_k^-1 are circulant kxk matrices (because R_k is one, and
>> R_k is positive definite => invertible)
> 
> Unfortunately they are not. Altough I'm willing to do that approximation
> (because R is "almost circulant"), if it considerably helps in this
> situation.

I made an error here. R is _not_ circulant and it isn't even close. Instead
it _is_ toeplitz. Sorry about that.

-- 
----
Jani Huhtanen
Tampere University of Technology, Pori

Jani,

> And to be more "precise": R is approximately toeplitz.

This is getting worse by the minute! I think I'm bailing out ... :-)

Best regards,
Andor

Jani Huhtanen wrote:

> Jani Huhtanen wrote:
>> abariska@student.ethz.ch wrote:
>> 
>>> ...
>>> 1. R_k and R_k^-1 are circulant kxk matrices (because R_k is one, and
>>> R_k is positive definite => invertible)
>> 
>> Unfortunately they are not. Altough I'm willing to do that approximation
>> (because R is "almost circulant"), if it considerably helps in this
>> situation.
> 
> I made an error here. R is _not_ circulant and it isn't even close.
> Instead it _is_ toeplitz. Sorry about that.
> 

And to be more "precise": R is approximately toeplitz.
-- 
----
Jani Huhtanen
Tampere University of Technology, Pori

> I know how to calculate the e because both R and A*R*A^T are invertible and
> there's no unknowns (except e). This is not the problem.
>
> The problem is that the calculation of e is a way too expensive operation
> and I need a another way to do it faster, even with the expense of
> accurary. Thus some approximation of e would be great.

Sorry, I cant help you Jani. Although I have worked with these
equations before, r was always unknown in my case, and you could use
the power method to compute both r and e...

http://college.hmco.com/mathematics/larson/elementary_linear/4e/shared/downloads/c10s3.pdf

I imagine in your case that if you work through a simple 3x3 example by
hand, you will quickly see if there is any structure to be exploited.

Porterboy

Hi Jani,

It can be shown that the expression is unity rank. I have slightly
reduced the problem to a simpler one though you would still require an
inverse of a matrix.

http://www.geocities.com/weech1at/repository/051111.pdf

Alan T

Alan Tan wrote:

> Hi Jani,
> 
> It can be shown that the expression is unity rank. I have slightly
> reduced the problem to a simpler one though you would still require an
> inverse of a matrix.
> 
> http://www.geocities.com/weech1at/repository/051111.pdf
> 
> Alan T

This is great! Thanks. 

a note to self: I really have to develop a good intuition for dealing with
submatrices.

-- 
----
Jani Huhtanen
Tampere University of Technology, Pori

Hi all!

Alan Tan showed how my equation can be simplified. Big thanks for him! 

Particularly intresting is the equation for R^-1. It gives R_{k}^-1 when
R_{k-1}^-1 is known. If this is applied recursively starting from say 1x1
matrix (scalar) one can get inverse of R with O(k^3) complexity (with 2*k^3
+ 3*k^2 + k mults), where k is the order (or size of the matrix R). 

Isn't this of the same complexity as inversion of symmetric positive
definite matrix through Cholesky decomp? 

If so, then it seems that I'm able to determine the "correct" order for the
linear predictor for "free" because I can check the minimum error at every
step of the algorithm.

Probably I have overlooked something or there is some other method for order
determination which is faster than this. Please tell me if you know such
method.

-- 
Jani Huhtanen
Tampere University of Technology, Pori

Jani Huhtanen wrote:
>...
>If this is applied recursively starting from say 1x1
> matrix (scalar) one can get inverse of R with O(k^3) complexity (with
> 2*k^3 + 3*k^2 + k mults), where k is the order (or size of the matrix R).

This is wrong. If I calculated correctly the amount of mults is given by
sum(k^2) from k=2 to K, where K is the final order (size of R). So the
correct amount of multiplications is (k^3)/3 + (k^2)/2 + k/6 - 1.

> 
> Isn't this of the same complexity as inversion of symmetric positive
> definite matrix through Cholesky decomp?

According to Golubs and van Loans "Matrix Computations 2nd ed." Cholesky
decomposition is O(n^3) and requires n^3/3 flops. This n^3/3 is an
approximation and it contains also additions. 

As a conclusion I think that the performance of the recursive calculation of
the R^-1 is comparable to Chol.decomp. Also if Alan Tan didn't do any
mistakes in his proof, this new (?) algorithm works for all symmetric
matrices (Am I right Alan?).  

Now the question is: Is this a _new_ algorithm? 
And is there even better ways to accomplish this (ie. finding the mse
associated with each predictor order and finding the R^-1 of the "optimal"
order)?

-- 
Jani Huhtanen
Tampere University of Technology, Pori

Hi Jani,

I should be quite correct in the formulation although personally I do
not think that this is a particularly new algorithm. For more
information about the inverse of the block matrices, you can refer to
the wonderful compilation of matrix operations by Roweis at
http://www.cs.toronto.edu/~roweis/notes/matrixid.pdf

Alan T


Jani Huhtanen wrote:
> Jani Huhtanen wrote:
> >...
> >If this is applied recursively starting from say 1x1
> > matrix (scalar) one can get inverse of R with O(k^3) complexity (with
> > 2*k^3 + 3*k^2 + k mults), where k is the order (or size of the matrix R).
>
> This is wrong. If I calculated correctly the amount of mults is given by
> sum(k^2) from k=2 to K, where K is the final order (size of R). So the
> correct amount of multiplications is (k^3)/3 + (k^2)/2 + k/6 - 1.
>
> >
> > Isn't this of the same complexity as inversion of symmetric positive
> > definite matrix through Cholesky decomp?
>
> According to Golubs and van Loans "Matrix Computations 2nd ed." Cholesky
> decomposition is O(n^3) and requires n^3/3 flops. This n^3/3 is an
> approximation and it contains also additions.
>
> As a conclusion I think that the performance of the recursive calculation of
> the R^-1 is comparable to Chol.decomp. Also if Alan Tan didn't do any
> mistakes in his proof, this new (?) algorithm works for all symmetric
> matrices (Am I right Alan?).
>
> Now the question is: Is this a _new_ algorithm?
> And is there even better ways to accomplish this (ie. finding the mse
> associated with each predictor order and finding the R^-1 of the "optimal"
> order)?
> 
> -- 
> Jani Huhtanen
> Tampere University of Technology, Pori