[SciPy-User] circulant matrices
Mon Feb 1 16:48:58 CST 2010
On Jan 31, 11:53 pm, Warren Weckesser <warren.weckes...@enthought.com>
> 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.
> 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
Thanks for the reference and the code again.
> > Any pointers would be appreciated.
> > Thanks,
> > Pål
> > _______________________________________________
> > SciPy-User mailing list
> > SciPy-U...@scipy.org
> [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)
More information about the SciPy-User