Linear operators for Julia

[Dominique Orban]

I just put together a basic linear operator package for Julia: https://github.com/dpo/linop.jl

The central idea is that operators act like matrices, but it would be too expensive to combine them before they are actually applied to a vector. Another benefit is that you can manipulate them as if they were matrices and write expressions such as y = (A' * A + B) * x without worrying that A' * A + B will actually be formed (it won't).

There are basic tests, but more should be added, and lots is still missing (e.g., block operators, parallel operators, ...). Entire algorithms could be wrapped as operators so that conceptually, you could write things like: x = LSQR(A) * b to apply the LSQR algorithm to solve the least-squares problem with A and b.

Feedback most welcome.

Cheers.

[Pierre-Yves Gérardy]

Rather than accumulating callbacks, couldn't you keep an array of operators and devectorize the computation? 

A potential drawback would be a loss of extensibility...

You could also cache the realized operator after the first evaluation (memoization).

—Pierre-Yves

[Pierre-Yves Gérardy]

Now that I think of it, wouldn't it make sense to have this mechanism for matrix operations in general (for Julia, that is)? 

You'd need immutable matrices. Only evaluate the operations when a result is indexed, using some kind of linear operation interpreter works iteratively without allocating temporary results where it makes sense (addition, substraction, transposition, scalar operations, maybe others)?

proposed syntax:

immutable_matrix = [| 1, 2 |

                    | 3, 4 |]

inline = [| 1, 2 | | 3, 4 |]

The use of the pipe is not ambiguous (bracket-pipe and value-pipe-space-pipe-value are not valid at the moment), and, by using comas, you don't have the pitfalls of significant white space as with hcat.

They could coexist for a while with their mutable cousins before they're deprecated.

I've seen the immutable Transpose proposition, but if the underlying 

Also:

threedee = [|| 1, 2 | | 5, 6 ||

            || 3, 4 | | 7, 8 ||]

fourdee = [||| 1, 2 | | 5, 6 ||

            || 3, 4 | | 7, 8 |||

           ||| a, b | | e, f ||

            || c, d | | g, h |||]

:-)

[Pierre-Yves Gérardy]

I forgot to finish this:

I've seen the immutable Transpose proposition, but if the underlying ...

 ... array is mutable, it is not very robust.

 —Pierre-Yves