ENH: linalg: add the function 'solve_circulant' for solving a circulant system

[WarrenWeckesser]

No description provided.

[coveralls]

Coverage remained the same when pulling 94cc658 on WarrenWeckesser:enh-linalg-solve-circ into 006a95f on scipy:master.

[jaimefrio]

It may require thinking about the UI a little harder, but I don't think you need to enforce 1-dimensionality. Instead, you could accept cs of shape (..., M) and bs of shape (..., M) or (..., M, N). If you can identify where in c and b is the dimension of shape M, you would simply have to take the fft along that axis, and let broadcasting do its magic. That way you could solve not only for N different RHS, but also for P different LHS by passing c of shape (P, M) and b of (M, N)... Sort of what happens in np.linalg.solve in 1.8.

[WarrenWeckesser]

Definitely something like that is possible. I thought about it while adding that strict dimension check, but decided to keep it simple to start. I've only spent a little time with the generalized API of solve, and it feels a bit awkward. Wouldn't it have made more sense to generalize the shape of b to (M, ...), which then includes (M,) and (M, N)? (Just roll M to the end for the purposes of broadcasting.) Anyway, I'll look at this some more, and see what makes sense for solve_circ.

[pv]

Generalized ufuncs roll from axes from the end, and moreover it's more efficient for C-order arrays.
Anyway, scipy.linalg should in any case follow the convention in numpy.linalg.

[jaimefrio]

Well, what would make sense from a gufunc broadcasting point of view is to have a single function that solved things like c = (.., M) and b = (..., M), and if you wanted to solve for multiple RHS pass in something like c = (M) and b = (N, M), if you wanted multiple LHS pass in c = (P, M) and b = (M), and if you wanted both then go with c = (P, 1, M) and b = (N, M), which would give you a return x = (P, N, M) where item x[p, n] would be the solution with the p-th LHS and the n-th RHS. But GESV follows the opposite convention for the RHS, and so the efficient handling of multiple RHS requires that b = (M, N), and that's what solve has been doing for ages, and sticking to the convention basically ruins having a nice consistent interface. If you look at the Python source code of np.linalg.solve, there are two different gufuncs called from solve, one for single RHS (solve1), other for multiple(solve), and the logic to select one or the other is not entirely satisfactory.

[argriffing]

Looks good to me. Maybe add tests with negative/imaginary/random c. Also it wouldn't bother me too much to type out solve_circulant.

[WarrenWeckesser]

@argriffing: Some more test cases is a good idea--I'll add some. I could live with the longer name if most folks prefer it.

[pv]

I think the longer name is justified, since then we'd have circulant and solve_circulant.

[WarrenWeckesser]

I've update the PR with more tests, and solve_circ is now solve_circulant.

To help with the discussion of the API, I also added the function solve_circulant_b. If there is agreement on the API, this will replace solve_circulant.

In the new API, there is no change to the basic case, where c and b are 1-d. For the case where either c or b is multidimensional, the arguments caxis and baxis (resp.) control the axis of the array that holds the coefficients. For broadcasting the computation, these axes are "ignored". For example, suppose c has shape (4,8,1) and caxis=1. We interpret this as a collection of length 8 vectors. The collection is an array with shape (4, 1). Simarly, an array b with shape (8, 3) and baxis=0 is a collection of length 8 vectors; the collection is an array with shape (3,). For broadcasting, the axes that hold the coefficients (which must be the same length) are ignored; the shapes of the collections must be compatible. In this example, the shapes of the collections are (4, 1) and (3,), which are compatible and produce a result with shape (4,3). So the result of solve_circulant_b will be a collection with shape (4,3) of length 8 vectors. With the default outaxis=-1, the result array has shape (4, 3, 8). (With outaxis=0, the result would have shape (8, 4, 3).)

This idea is quite general, and could be used elsewhere. It generalizes the API of the gufuncs in numpy.

Note that the special case where b has shape (M, N) is dropped in solve_circulant_b. In that case, you would use baxis=0.

[coveralls]

Coverage remained the same when pulling 64f0daf on WarrenWeckesser:enh-linalg-solve-circ into 006a95f on scipy:master.

[ev-br]

Bump.
Would be good to revive this in time for scipy 0.15, I'd think.

[pv]

The _b version API is in principle fine from broadcasting POV, but I'm a
bit worried about the surprise in solve_circulant(c, b) where b is a
matrix (solves transposed problem compare to solve()). Documentation can
of course reduce the surprise, but on the other hand picking the ufunc
based on the number of dimensions as in solve() is maybe not that bad
either?

[argriffing]

I haven't read all of the API discussion on this PR but I'd expect solve_circulant(c, b) to be a specialized implementation of solve(circulant(c), b).

[WarrenWeckesser]

I updated the pull request. The proposed signature of solve_circulant is now

def solve_circulant(c, b, singular='raise', tol=None, caxis=-1, baxis=0, outaxis=0):

When c is a 1-d array, @argriffing's suggested behavior is followed: solve_circulant(c, b) is a specialized implementation of solve(circulant(c), b), where solve is scipy.linalg.solve.

Note that scipy.linalg.solve and numpy.linalg.solve do not handle high dimensional b arrays the same way. If, for example, a has shape (3, 3) and b has shape (3, 2, 5), scipy.linalg.solve(a, b) returns an array with shape (3, 2, 5); numpy.linalg.solve(a, b) raises a ValueError. solve_circulant follows the convention of scipy.linalg.solve.

[rgommers]

This looks ready. Besides a rebase, is there anything left to do here?

[WarrenWeckesser]

Rebased. I don't have any more changes in mind.

[rgommers]

@pv @argriffing maybe time for a final review and then green button it?

[argriffing]

Yes, it has tests asserting that solve_circulant(c, b) is almost equal to solve(circulant(c), b), including tests for the 2d rhs mentioned by @pv, so I think it is safe to merge.

[argriffing]

Thanks @WarrenWeckesser sorry it took a year to merge :/