toeplitz matrix [Matlab]

Prev: LaTeX interpreter problem in R2009a?
Next: S-function input data type as enumeration

From: Xiaoyi on 11 May 2010 05:27

"Matt J " <mattjacREMOVE(a)THISieee.spam> wrote in message <hlivn3$5he$1(a)fred.mathworks.com>...
> "Bc " <bcortez23(a)gmail.com> wrote in message <hlikch$6o5$1(a)fred.mathworks.com>...
> > Hi all,
> >
> > I am trying to create a toeplitz matrix which will be of a really large size. In fact, when doing a symmetric for r = 9216x1
> >
> > toeplitz(r)
> >
> > I get an error. The worst part is that this is a small r for what I need, so I wanted to ask what there is in terms of a solving Ax = b for x, given b and that A is a toeplitz matrix. I was going to do inv(A) or A\b, but I can't even create the toeplitz in the first place.
> ===============
>
> A little more info is needed on how sparse your matrix A is. If it is sparse, my interpMatrix tool might be of help:
>
> http://www.mathworks.com/matlabcentral/fileexchange/26292-regular-control-point-interpolation-matrix-with-boundary-conditions
>
> It is fairly well optimized at creating sparse Toeplitz and Toeplitz-like matrices.
>
> If A is not sparse, it is no surprise that it is exceeding memory limits. In that case, you can probably exploit the fact that A*x=conv(A(:,1),x) . If you zero-pad this linear convolution so as to make it a circulant convolution, your equation has the closed form solution
>
> fft(x)=fft(b)./fft(A(:,1))
>
> assuming that fft(A(:,1)) is non-zero everywhere. Otherwise, there is no solution.

Can you please elaborate on this part? I tried it and it did not work. Maybe I was doing something wrong? Suppose A is a symmetric Toeplitz matrix (dense), how would you solve Ax=b with fft and ifft?

Thanks!

From: Matt J on 11 May 2010 11:12

"Xiaoyi " <shawREMOVE_liu(a)hotmail.THIScom> wrote in message <hsb7t7$pls$1(a)fred.mathworks.com>...

> > If A is not sparse, it is no surprise that it is exceeding memory limits. In that case, you can probably exploit the fact that A*x=conv(A(:,1),x) . If you zero-pad this linear convolution so as to make it a circulant convolution, your equation has the closed form solution
> >
> > fft(x)=fft(b)./fft(A(:,1))
> >
> > assuming that fft(A(:,1)) is non-zero everywhere. Otherwise, there is no solution.
>
> Can you please elaborate on this part? I tried it and it did not work. Maybe I was doing something wrong? Suppose A is a symmetric Toeplitz matrix (dense), how would you solve Ax=b with fft and ifft?
========

If A is Toeplitz, there is some vector, k, such that

A*x=conv(k,x,'same')

The equation you must solve is therefore

conv(k,x,'same')=b

If k is zero padded to the same length as b, then this can also be represented as a cyclic convolution and taking FFTs gives

FFT(k)*FFT(x)=FFT(b)

leading to

x=IFFT(FFT(b)./FFT(k))

From: Xiaoyi on 11 May 2010 18:44

"Matt J " <mattjacREMOVE(a)THISieee.spam> wrote in message <hsbs45$q20$1(a)fred.mathworks.com>...
> "Xiaoyi " <shawREMOVE_liu(a)hotmail.THIScom> wrote in message <hsb7t7$pls$1(a)fred.mathworks.com>...
>
> > > If A is not sparse, it is no surprise that it is exceeding memory limits. In that case, you can probably exploit the fact that A*x=conv(A(:,1),x) . If you zero-pad this linear convolution so as to make it a circulant convolution, your equation has the closed form solution
> > >
> > > fft(x)=fft(b)./fft(A(:,1))
> > >
> > > assuming that fft(A(:,1)) is non-zero everywhere. Otherwise, there is no solution.
> >
> > Can you please elaborate on this part? I tried it and it did not work. Maybe I was doing something wrong? Suppose A is a symmetric Toeplitz matrix (dense), how would you solve Ax=b with fft and ifft?
> ========
>
> If A is Toeplitz, there is some vector, k, such that
>
> A*x=conv(k,x,'same')
>
> The equation you must solve is therefore
>
> conv(k,x,'same')=b
>
> If k is zero padded to the same length as b, then this can also be represented as a cyclic convolution and taking FFTs gives
>
> FFT(k)*FFT(x)=FFT(b)
>
> leading to
>
> x=IFFT(FFT(b)./FFT(k))

I am still confused. For a simple example, let's say
A=toeplitz([5 2 1]);
b=[1;1;1];
how sould you find such a vector k (3x1) and how would you zero-pad it to make it circular and then solve the system Ax=b?

From: Xiaoyi on 11 May 2010 18:46

From: Matt J on 12 May 2010 14:51

"Xiaoyi " <shawREMOVE_liu(a)hotmail.THIScom> wrote in message <hscmnf$f4r$1(a)fred.mathworks.com>...

>
> I am still confused. For a simple example, let's say
> A=toeplitz([5 2 1]);
> b=[1;1;1];
> how sould you find such a vector k (3x1) and how would you zero-pad it to make it circular and then solve the system Ax=b?
==============

Sorry, I don't think it will work after all, because b only contains a truncation of the convolution represented by A*x and not the whole convolution.

However, here is an expanded discussion of what I was getting at, using your example above. First, the solution to this equation is

>> x=A\b

x =

0.1364
0.0909
0.1364

Note that if I define the vector k=[1 2 5 2 1], then the equation is equivalent to
conv(x,k,'same')=b

>> conv(x,k,'same') % = A*x = b

ans =

1
1
1

I would like to invert this equation using FFTs, but because of the truncation induced by the 'same' option, I cannot. If, however, I had the more complete data

>> b=conv(x,k)

b =

0.1364 0.3636 1.0000 1.0000 1.0000 0.3636 0.1364

then I could invert the convolution using FFTs as follows

>> N=length(b); xx=ifft(fft(b,N)./fft(k,N),N); xx=xx(1:3).'

xx =

0.1364
0.0909
0.1364

|
Pages: 1
Prev: LaTeX interpreter problem in R2009a?
Next: S-function input data type as enumeration