What does inv(A) do when A is a Toeplitz Matrix?
Elvis Chen
greetings,
I'm trying to convert some of my matlab codes into C++ equivalent program.
For matrix inversion, I used LU decomposition with back-substitution to
find the inverse of a give matrix. However, when the matrix is large
(> 2000x2000) and is a Toeplitz matrix, it is very slow.
However, Matlab can find the inverse of such matrices without any
difficulty. On an Ultra-60 with 2gig RAM, it finds the inverse in just
about 10 seconds.
My question is, how does matlab do it? Is there a special algorithm to
find the inverse of a toeplitz matrix?
thanx in advance,
Elvis
Jonathan Shekter
>difficulty. On an Ultra-60 with 2gig RAM, it finds the inverse in just
>about 10 seconds.
>
>My question is, how does matlab do it? Is there a special algorithm to
>find the inverse of a toeplitz matrix?
Not that I know of. If such an algorithm existed, it would be of tremendous
value to the signal processing community. Many filtering algorithms involve
large Toeplitz matrices, and such matrices are enough of a problem that at the
moment the standard technique is to instead approximate this matrix by a
circulant matrix. Since all circulant matrices are diagonalized by a FFT,
inverses and many other operations can be computed quickly on circulant
matrices.
Aaanyway, I don't know how MATLAB does it but I have theories:
1) Often, Toeplitz matrices are sparse and in this case they are usually also
band-diagonal. There exist very fast (linear time!) algorithms for solution of
band-diagonal systems.
2) Since MATLAB is a mature product whose primary function is the manipulation
of matrices, it may just be that it uses standard algorithms but it has
bitchin' matrix inversion routines which have been optimized within an inch of
their lives.
3) Getting way out there: possibly since circulants are ''close'' to Toeplitz
matrices, one could use an iterative algorithm which starts with the inverse of
the corresponding circulant matrix as an initial guess.
Bergers
>From: Elvis Chen ch...@cs.queensu.ca
>Date: 4/1/00 6:24 PM Mountain Daylight Time
>Message-id: <Pine.SOL.3.95.1000401192004.17864B-100000@apex>
>
>
>greetings,
>
>I'm trying to convert some of my matlab codes into C++ equivalent program.
>For matrix inversion, I used LU decomposition with back-substitution to
>find the inverse of a give matrix. However, when the matrix is large
>(> 2000x2000) and is a Toeplitz matrix, it is very slow.
>
>difficulty. On an Ultra-60 with 2gig RAM, it finds the inverse in just
>about 10 seconds.
>
>My question is, how does matlab do it? Is there a special algorithm to
>find the inverse of a toeplitz matrix?
>
>
>Elvis
You might want to look into the Levinson (-Durbin) or Schur algorithms for
solving a system of equations involving a Toeplitz matrix. These Toeplitz
specific algorithms only require on the order of n^2 operations versus n^3.
Scott
News Groups
solving
toeplitz(1:512)\(1:512)'
takes my MATLAB 4.6 seconds, while solving it with our (compiled) levinson
routine, takes
10 milli-seconds!
MATLAB has the LEVINSON (-DURBIN) routine, which can be used to solve
toeplitz equations if the
RHS is of a certain form. But you will need to write your own generalised
LEVINSON routine.
Brett.
"Elvis Chen" <ch...@cs.queensu.ca> wrote in message
news:Pine.SOL.3.95.1000401192004.17864B-100000@apex...
Elvis Chen
> toeplitz(1:512)\(1:512)'
> takes my MATLAB 4.6 seconds, while solving it with our (compiled) levinson
> routine, takes
> 10 milli-seconds!
> > I'm trying to convert some of my matlab codes into C++ equivalent program.
> > For matrix inversion, I used LU decomposition with back-substitution to
> > find the inverse of a give matrix. However, when the matrix is large
> > (> 2000x2000) and is a Toeplitz matrix, it is very slow.
thanx,
perhaps I'm missing something. Over the past few days, I have coded
Gauss-Jordan Elimination and LU-decomposition to find the inverse of a
matrix in C++. The code were originally taken from Numerical Recipes,
which I assumed was pretty good and efficient. However, when I tried to
use them to find the inverse of a matrix with size 2048x2048, even on an
Ultra 60 with 2gig RAM it takes forever! On the same machine, matlab 5.2
only takes 10 seconds! If what you said is true, that Matlab doesn't care
if a matrix is Toeplitz, then perhaps using Numerical Recipes was a wrong,
very wong, starting point in the first place.
I'll look into Levison-Durbin algorithm. Does anyone know if I can find a
free code of it?
thanx again,
Elvis
Johan Kullstam
> On Mon, 3 Apr 2000, News Groups wrote:
>
> > toeplitz(1:512)\(1:512)'
> > takes my MATLAB 4.6 seconds, while solving it with our (compiled) levinson
> > routine, takes
> > 10 milli-seconds!
> > > I'm trying to convert some of my matlab codes into C++ equivalent program.
> > > For matrix inversion, I used LU decomposition with back-substitution to
> > > find the inverse of a give matrix. However, when the matrix is large
> > > (> 2000x2000) and is a Toeplitz matrix, it is very slow.
>
> thanx,
>
> perhaps I'm missing something. Over the past few days, I have coded
> Gauss-Jordan Elimination and LU-decomposition to find the inverse of a
> matrix in C++. The code were originally taken from Numerical Recipes,
> which I assumed was pretty good and efficient.
numerical recipies is generally rather weak. their fortran isn't half
bad, but their C code is truly awful.
> However, when I tried to
> use them to find the inverse of a matrix with size 2048x2048, even on an
> Ultra 60 with 2gig RAM it takes forever! On the same machine, matlab 5.2
> only takes 10 seconds! If what you said is true, that Matlab doesn't care
> if a matrix is Toeplitz, then perhaps using Numerical Recipes was a wrong,
> very wong, starting point in the first place.
>
> I'll look into Levison-Durbin algorithm. Does anyone know if I can find a
> free code of it?
you might look at the lattice filter stuff. the lattice filter is
based on L-D.
--
J o h a n K u l l s t a m
[kull...@ne.mediaone.net]
Don't Fear the Penguin!
Tom Krauss
Here is a c-mex source routine to solve a toeplitz system
efficiently. I got it from one of my favorite DSP books,
Roberts and Mullis. Drawback: it only works for real data.
If you extend it to the complex case please share it with the
rest of us.
Cheers,
Tom Krauss
PhD Student
--------------------- m-file help function:
toepsolve.m ---------------------
%function h = toepsolve(r,q);
%TOEPSOLVE Solve Toeplitz system of equations.
% Solves R*h = q, where R is the symmetric Toeplitz matrix
% whos first column is r
% Assumes all inputs are real
% Inputs:
% r - first column of Toeplitz matrix, length n
% q - rhs vector, length n
% Outputs:
% h - length n solution
%
% Algorithm from Roberts & Mullis, p.233
%
% Author: T. Krauss, Sept 10, 1997
--------------------- c-mex source: toepsolve.c ------------------------
#include <math.h>
#include "mex.h"
/* function h = toepsolve(r,q);
* TOEPSOLVE Solve Toeplitz system of equations.
* Solves R*h = q, where R is the symmetric Toeplitz matrix
* whos first column is r
* Assumes all inputs are real
* Inputs:
* r - first column of Toeplitz matrix, length n
* q - rhs vector, length n
* Outputs:
* h - length n solution
*
* Algorithm from Roberts & Mullis, p.233
*
* Author: T. Krauss, Sept 10, 1997
*/
void mexFunction(
int nlhs,
mxArray *plhs[],
int nrhs,
const mxArray *prhs[]
)
{
double *a,*h,beta;
int j,k;
double eps = mxGetEps();
int n = (mxGetN(prhs[0])>=mxGetM(prhs[0])) ? mxGetN(prhs[0]) :
mxGetM(prhs[0]) ;
double *r = mxGetPr(prhs[0]);
double *q = mxGetPr(prhs[1]);
double alpha = r[0];
n = n - 1;
plhs[0] = mxCreateDoubleMatrix(n+1,1,0);
h = mxGetPr(plhs[0]);
h[0] = q[0]/r[0];
a = mxCalloc((n+1)*(n+1),sizeof(double));
if (a == NULL) {
mexErrMsgTxt("Sorry, failed to allocate buffer.");
}
a[(0*(n+1))+0] = 1.0;
for (k = 1; k <= n; k++) {
a[(k*(n+1))+k-1] = 0;
a[(0*(n+1))+k] = 1.0;
beta = 0.0;
for (j = 0; j <= k-1; j++) {
beta += r[k-j]*a[(j*(n+1))+k-1];
}
beta /= alpha;
for (j = 1; j <= k; j++) {
a[(j*(n+1))+k] = a[(j*(n+1))+k-1] - beta*a[((k-j)*(n+1))+k-1];
}
alpha *= (1 - beta*beta);
h[k] = q[k];
for (j = 0; j <= k-1; j++) {
h[k] -= r[k-j]*h[j];
}
h[k] /= alpha;
for (j = 0; j <= k-1; j++) {
h[j] += a[((k-j)*(n+1))+k]*h[k];
}
} /* loop over k */
mxFree(a);
return;
}
Elvis Chen
> Here is a c-mex source routine to solve a toeplitz system
> efficiently. I got it from one of my favorite DSP books,
> Roberts and Mullis. Drawback: it only works for real data.
> If you extend it to the complex case please share it with the
> rest of us.
hi Tom,
I was given a piece of Trench algorithm that computes the INVERSE of a
Toeplitz matrix directly. There are two versions: one works for
symmetric Toeplitz matrices, the other is a generalized version that works
for non-symmetric toeplitz matrix.
However, the generalized version seems to give the wrong answer when the
given Toeplitz matrix is NON-symmetrical. It does give the right answer
then the given Toeplitz matrix is symmetrical.
I have since translated them into C++ program (not cmex), and I'm willing
to share it.
The problem is, the person who gave me the printed copy lost it's
citation, so I don't know if such code is from a book, someone's paper, or
someone's thesis. That means, I don't know for certain it someone is
holding the copy-right of it. Assuming it isn't violating any law, I'm
willing to send anyone a copy if I'm contacted by private email.
ttyl,
Elvis
Anitha Florence Vinola
> >However, Matlab can find the inverse of [large Toeplitz] matrices without any
> >difficulty. On an Ultra-60 with 2gig RAM, it finds the inverse in just
> >about 10 seconds.
> >
Greg Heath
wrote:
> jshek...@interlog.com (Jonathan Shekter) wrote in message <J2HF4.110044$Hq3.2787...@news2.rdc1.on.home.com>...
>
> - Show quoted text -
inv(toeplitz(n)) = 1/n
Hope this helps.
Greg