[SciPy-User] circulant matrices

[Paul]

On Jan 31, 11:53 pm, Warren Weckesser <warren.weckes... at enthought.com>
wrote:
> Paul wrote:
> > I would like to find,
>
> > argmin_c   norm( x - dot(phi, c), 2)
>
> > where x, c are vectors and phi is a circulant matrix. Is there a way
> > to do this with a built-in numpy or scipy function? LAPACK doesn't
> > seem to have a circulant matrix type. There is a Fortran library
> > NAPACK on netlib.org that has a circulant matrix type though I don't
> > know anything about it. I guess I could try to add these functions
> > with f2py.
>
> > I tried to solve this problem using a naive deconvolution approach
> > with fft's but it's not robust. I guess I have to do some more
> > research. I can't really afford to treat phi as anything other than
> > circulant (like a general banded or dense matrix) as it has to be
> > fast.
>
> Is phi singular?  If not, then wouldn't you just solve dot(phi,c) = x for c?

Thank you for your example code! It helped me realize that I have been
doing something wrong.

phi is not singular and you are right that you can solve dot(phi,c) =
x directly, with the fast implementation that you provided. I am
guessing phi not being singular also means fc = fft(c) has no zeros.

My mistake is that my phi is not exactly circulant, that is, it has
several additional rows. That is why I thought my problem is
overdetermined and required an optimization. I am representing a
correlation of a vector phi with c (where len(phi) << len(c) and I
need to correct my padding so that I can solve it in the way you
suggested. In the full problem, the correlation is actually 2-d,
complicating things a little further.

>
> Based on the wikipedia articlehttp://en.wikipedia.org/wiki/Circulant_matrix,
> I implemented a quick test of the FFT method--see the attached script.  
> (Sorry,
> my notation is the usual Ax=b instead of phi*c=x.)  It seems to work, and it
> is much faster than linalg.solve(), but the script is the full extent of my
> experience with circulant matrices (i.e. probably even more naive than
> what you have already tried), so there are probably details that I am
> missing.

Thanks for the reference and the code again.

>
> Warren
>
> > Any pointers would be appreciated.
>
> > Thanks,
> > Pål
> > _______________________________________________
> > SciPy-User mailing list
> > SciPy-U... at scipy.org
> >http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
>
> [circtest.py< 1K ]import time
> import numpy as np
> from scipy import fft, ifft
>
> # Solve A*x = b when A is circulant.
>
> # Create a random circulant matrix.
> n = 512
> c = 2*np.random.random(n) - 1.0
> A = np.empty((n,n), dtype=float)
> for k in range(n):
>     A[:,k] = np.roll(c,k)
>
> # Create the b vector.
> b = np.arange(n)
>
> # First use linalg.solve() for comparison.
> t = time.time()
> x1 = np.linalg.solve(A,b)
> dt1 = time.time() - t
> print "Using solve: %.3e seconds" % dt1
>
> # Next use the FFT.
> # Note that we only need c here; A was created to be used in solve().
> t = time.time()
> fb = fft(b)
> fc = fft(c)
> x2 = ifft(fb / fc)
> dt2 = time.time() - t
> print "Using FFT:   %.3e seconds" % dt2
>
> print "Speedup: %.3g" % (dt1/dt2)
>
> # Verify that the solutions are the same.
> assert np.allclose(x1, x2.real)
>

Pål.