Change fft(x) default behavior for matrices

[stevengj]

Somewhat unexpectedly (especially for people used to Matlab), Julia gives:

julia> fft([1. 2. 3.])
1x3 Complex128 Array:
 1.0+0.0im  2.0+0.0im  3.0+0.0im

because Julia defines fft(X) = fft(X, 1), whereas what the user wants is probably

julia> fft([1. 2. 3.], 2)
1x3 Complex128 Array:
 6.0+0.0im  -1.5+0.866025im  -1.5-0.866025im

You might consider making the fft (and ifft and bfft and rfft) functions transform the second dimension when they act on a row vector. (The user can always do fft(x,1) if they want to force the current behavior.)

Update: See below for better suggestions.

[andreasnoack]

After getting used to Julian vectors I like the present behaviour. If you define a vector instead of a matrix you get the expected

julia> fft([1.,2.,3.])
3-element Complex128 Array:
       6.0+0.0im
 -1.5+0.866025im
 -1.5-0.866025im

[stevengj]

Good point, a lot of it is just that I'm not used to commas.

[timholy]

This has been discussed before in previous contexts, but I can't find the links now, so... while I came from Matlab, I'd declare Matlab's behavior here to be completely broken (I felt that way long before discovering Julia). Say I write a routine that expects a matrix input, and I use fft intending it to compute along columns. Now, one time out of 10^5, some routine feeds it a matrix with just one row. Suddenly my routine switches to computing the fft along rows. Big, nasty, difficult-to-detect bug.

Personally I vote we make it an error to call fft on anything but a vector without the dimension index.

[stevengj]

I agree with @timholy ... defaulting fft(x) = fft(x,1) for matrices is just asking for confusion.

If you want to break with Matlab's confusing conventions, I'd argue for fft(x) = fftn(x) by default: i.e. it does a multidimensional FFT unless you specify a dimension index explicitly fft(x,3). And for that matter I would support multiple dimension indices, so that e.g. fft(x,[1,3]) transformed dimensions 1 and 3 of a multidimensional array. And for that matter I wouldn't bother exporting fftn, fft2, etc.

[timholy]

I thought about fft = fftn by default too. I'm just scared about the consequences for people who port over Matlab code. I guess my philosophy is "errors are OK" (they're easy to catch), but "incompatible behavior is more scary".

[stevengj]

Okay, that's reasonable to me.

Technically, in that case you still may have a problem, because Matlab's fft function defines fft(x,n) as "zero-pad to length n", whereas to transform dimension d you need to do fft(x,[],d). But arguably people doing fft(x,n) in Matlab are almost always going to have n > ndim(x), so Julia will give a bounds error.

[timholy]

Hmm. I hadn't realized that. I like the bounds error, but there's nothing really ideal.

If no one chimes in against this plan (specifically: deleting the default dimension for anything other than a vector input), I'll implement it by the end of the weekend.

[staticfloat]

I like it, makes sense and has basically no cost to the developer.

[JeffBezanson]

It would be great to eliminate the clutter of fft2, fft3, and fftn.

[stevengj]

Let me propose a slightly different alternative:

• fft(x) works only on 1d arrays, and throws an error for ndims(x) > 1
• fftn(x) works for multidim arrays, like Matlab
• ffts(dims, x) performs FFTs on dimension(s) dims (which can be an integer or an array of integers, and maybe a range too) of x [where "s" could stand for "subset" or the "s" in "FFTs", depending on your preference]

Defining a separate fft and fftn in this way is somewhat ugly, but avoids silently giving wrong results for code ported over from Matlab as per @timholy. Defining a separate ffts function in this way has several advantages:

• It eliminates any possibility of conflict with the Matlab fft(x,n) syntax.
• It works better with the proposed "optimized FFT" interface from issue #1805: fft(x) and fftn(x) have optimized variants plan_fft(x,flags) and plan_fftn(x,flags), while ffts(dims,x) has the optimized variant plan_ffts(dims,x,flags), and there is no possibility that dims will get mixed up with flags.

[ViralBShah]

I like the proposed interface. It does not confuse matlab users, and is more powerful. How about

ffts(x, dims)

since we are passing dims as the last argument to many functions in line with what matlab does. We already have the notion of a region with Dimspec for reductions as in sum(a, (2,3)), and the same can be carried over to fft.

[stevengj]

Okay. I'm playing with an implementation for this and issue #1805 , but probably won't have anything for a few days due to the holidays.

[stevengj]

Question: are there any circumstances in which an array may not be 16-byte aligned? e.g. is it possible to make an array from a subset of another array without a copy? Or from a non-aligned pointer returned by a C library?

And are you sure that arrays are 16-byte aligned? array.c seems to use jl_gc_managed_malloc, which just calls malloc. malloc is 16-byte aligned on MacOS X but is only 8-byte aligned on e.g. 32-bit x86 glibc.

[ViralBShah]

@JeffBezanson will have to say. I distinctly remember us making the change for DSFMT, which also requires 16-byte aligned memory for the same reasons.

[JeffBezanson]

I only use malloc above a certain size, where it seems the result is always aligned 16. I don't recall whether this is guaranteed, but in practice it seems to be. Subarrays are an issue. Functions that can work on subarrays and have strict alignment needs will probably have to check. When wapping pointers from c we could add a check as well. On Dec 22, 2012 9:46 AM, "Viral B. Shah" notifications@github.com wrote:

@JeffBezanson https://github.com/JeffBezanson will have to say. I distinctly remember us making the change for DSFMT, which also requires 16-byte aligned memory for the same reasons.

— Reply to this email directly or view it on GitHubhttps://github.com/JuliaLang/julia/issues/1806#issuecomment-11637724.

[stevengj]

@JeffBezanson Have you tested this in 32-bit mode? x86_64 malloc is guaranteed to be 16-byte aligned on most OS's, and in 32-bit mode on MacOS X , but only 8-byte alignment is guaranteed in 32-bit x86 on Win32 (and even then only for allocating at least 8 bytes: http://msdn.microsoft.com/en-us/library/83ythb65.aspx) or glibc (http://www.gnu.org/software/libc/manual/html_node/Aligned-Memory-Blocks.html#Aligned-Memory-Blocks).

Would be safer to use posix_memalign on Unices. See also the comments on various operating systems in FFTW's kernel/kalloc.c file ... you could also just use FFTW's fftw_malloc and fftw_free functions.

[andreasnoack]

In the proposed implementation fft could be overloaded instead of having all the variants with s and n. This is a general consideration that applies more widely that fft. Julia has several functions that behave different from Matlab, e.g. sum so wouldn't a general advise to read the documentation be better than adding function names just to be sure that people used to Matlab get an error?

[ViralBShah]

We can stick it all in fft once we have keyword arguments. Without that, it would be inconsistent with MATLAB usage, and generally confusing with so much different functionality crammed into one function.

[stevengj]

@ViralBShah I think @andreasnoakjensen was suggesting using fft(x) instead of fftn(x), and fft(x,dims) for ffts. This is similar to sum, and isn't so much functionality. Then you can get rid of fftn, fft2, fft3, and the proposed ffts.

I prefer this solution myself. The only question is whether one is worried about silently giving different results for Matlab users as per @timholy. But I am persuaded by the argument that you are already doing this for sum and similar functions. It would be more internally consistent for Julia to treat fft the same as sum, and once users learn of the different behavior for sum they will probably expect it for fft too anyway.

[johnmyleswhite]

I agree with @andreasnoackjensen and @stevengj that breaking from Matlab is less worrisome than breaking from Julia's own sum. I also believe that R's fft behaves according to their proposal.

[ViralBShah]

I too think that just using fft is the cleanest design.

sum without arguments gives a scalar, which is the sum over the entire array. So, matlab users will be surprised, but will discover right away what is going on, since most codes will break. In the case of fft, the code will keep working silently, which is the worry.

That said, in general we have taken the approach of having clean and meaningful interfaces as the first priority, and retaining similarity to matlab when it does not compromise the former. I am willing to go with the clean interface, and just seeing how this turns out.

The fact that it is similar to R's fft makes the argument more compelling.

[johnmyleswhite]

From the R docs for fft:

 When ‘z’ is a vector, the value computed and returned by ‘fft’ is
 the unnormalized univariate Fourier transform of the sequence of
 values in ‘z’.

 When ‘z’ contains an array, ‘fft’ computes and returns the
 multivariate (spatial) transform.  If ‘inverse’ is ‘TRUE’, the
 (unnormalized) inverse Fourier transform is returned, i.e., if ‘y
 <- fft(z)’, then ‘z’ is ‘fft(y, inverse = TRUE) / length(y)’.

 By contrast, ‘mvfft’ takes a real or complex matrix as argument,
 and returns a similar shaped matrix, but with each column replaced
 by its discrete Fourier transform.  This is useful for analyzing
 vector-valued series.

[timholy]

Just voicing the main negative of the change, but I agree there are good reasons in favor, as well. I'm fine with whatever those who actually do the work decide.

[StefanKarpinski]

For clarity of discussion can someone give an example of where a Matlab call to fft and the proposed Julia fft would not behave in the same way?

[ViralBShah]

This will behave differently:

a = rand(5,5)
b = fft(a)

In Matlab, that would compute the FFT of every column of a. In the new proposed interface, julia will instead compute the n-d FFT (fftn in Matlab) for an n-d array. In order to get the Matlab-like behaviour, you have to pass the optional dim argument (which behaves the same way as matlab).

[stevengj]

@StefanKarpinski: Suppose x is a matrix (a 2d array). The proposed fft(x) would compute what Matlab calls fftn(x), the 2d DFT of x.

Matlab's fft(x), in contrast, computes the 1d DFT of each column of x. In the proposed interface, this would be achieved in Julia by fft(x,2). [In Matlab, fft(x,2) would mean to zero-pad the columns of x to length 2, or truncate them if they have length > 2. This is mainly an artifact of the days when fft was slow for non power-of-two sizes.]

Matlab's interface makes some sense if you adopt the philosophy that fft(x) is equivalent F * x where F is the DFT matrix. For higher-dimensional matrices, you could interpret it as a tensor contraction with the first index. However, Matlab itself is inconsistent because fft(x) on a row-vector x computes the 1d DFTs of the row rather than the (trivial) 1d DFT of the columns (this is what Julia does now and is the source of this thread). [In my experience, many Matlab users accidentally employ fft(x) when what they want is fftn(x).]

Matlab also has a somewhat silly proliferation of interfaces: fft2, fft3, and fftn for 2d, 3d, and nd DFTs. These are all subsumed by my proposed interface. Note that fft2 is not equivalent to fftn: e.g. for a 3d array x, fft2(x) computes the 2d DFT of each slice of x, which is equivalent in my proposed interface to fft(x, (1,2)) (transform along the first two dimensions of x).

[StefanKarpinski]

Right. It sounds like making fft compute the n-dimensional fft would be best.

[stevengj]

Fixed by #1856.