# Unified Interface for Linear Solving

**URL:** <https://discourse.julialang.org/t/unified-interface-for-linear-solving/699>\
**Category:** Numerics\
**Created:** [December 2, 2016, 5:01pm UTC](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699 "2016-12-02T17:01:01Z")\
**Posts on this page:** 15\
**Page:** 2

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [December 11, 2016, 5:17am UTC](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699/21 "2016-12-11T05:17:27Z")

</div>

> [@barche](#):
>
> solve!(x, A, b, parameters)

In Base, this routine is called `A_ldiv_B!`

---

<div class="post-metadata">

**Author:** ![barche](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/barche/32/79_2.png) [@barche](https://discourse.julialang.org/u/barche)\
**Post date:** [December 11, 2016, 10:15am UTC](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699/22 "2016-12-11T10:15:11Z")

</div>

> [@ChrisRackauckas](#):
>
> The problem with solve is I don’t see how it can match this two-step approach in an easy manner and get the same efficiency. It can get lots of these things right, but the initialization is really crucial for many linear solving setups (not for iterative solvers, but just thinking about CPU-based iterative solvers is really close-minded here). And I don’t see why you’d want to avoid : MATLAB had some bad syntax choices, but \ is one choice that I find to be very clear and elegant.

You’re absolutely right that initialization is necessary, in Trilinos the “A” is also something called a “linear operator with solve”, which is typically an initialized factorization, so it would present the same workflow, except that there is the extra parameters option for passing parameters you might want to vary at each step and that don’t require recomputing the initialized operator, maybe things like the tolerance or the number of iterations.

That said, I like your example and it would fit my future use case (solving the Navier-Stokes equations) well, but I think it would be nice to also have a more advanced syntax that allows passing in extra parameters at each step if needed.

How standard is the `@into!` syntax? It is much nicer than having to remember a non-intuitive ordering of function arguments, I think it would be good to have macros like this in Base so they are more universally recognised.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 11, 2016, 6:04pm UTC](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699/23 "2016-12-11T18:04:52Z")

</div>

> [@barche](#):
>
> That said, I like your example and it would fit my future use case (solving the Navier-Stokes equations) well, but I think it would be nice to also have a more advanced syntax that allows passing in extra parameters at each step if needed.

We’d have to try to standardize some arguments if we want that to work nicely. Can you give an example of what you’re thinking of?

> [@barche](#):
>
> How standard is the @into! syntax? It is much nicer than having to remember a non-intuitive ordering of function arguments, I think it would be good to have macros like this in Base so they are more universally recognised.

`@into!` is just a nice way from InplaceOps.jl to use ` A_ldiv_B!`. ` A_ldiv_B!` is standardized. I prefer the \ look with a macro in front. Pretty sure this macro isn’t standard, but I would hope it would be become so.

---

<div class="post-metadata">

**Author:** ![barche](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/barche/32/79_2.png) [@barche](https://discourse.julialang.org/u/barche)\
**Post date:** [December 12, 2016, 8:02am UTC](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699/24 "2016-12-12T08:02:03Z")

</div>

> [@ChrisRackauckas](#):
>
> We’d have to try to standardize some arguments if we want that to work nicely. Can you give an example of what you’re thinking of?

I was thinking of just adding an extra argument `parameters` with an abstract type (`abstract SolveParameters` ?), so using multiple dispatch it could be specialized to contain anything.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 12, 2016, 8:08am UTC](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699/25 "2016-12-12T08:08:24Z")

</div>

> [@barche](#):
>
> I was thinking of just adding an extra argument parameters with an abstract type (abstract SolveParameters ?), so using multiple dispatch it could be specialized to contain anything.

What would that give us over the method using “God types” and closures? That’s the example I’m looking for.

---

<div class="post-metadata">

**Author:** ![barche](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/barche/32/79_2.png) [@barche](https://discourse.julialang.org/u/barche)\
**Post date:** [December 12, 2016, 10:26am UTC](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699/26 "2016-12-12T10:26:19Z")

</div>

You could set e.g. tolerance at each step:

```julia
function OneStep(f,...,linear_solver=lufact) 
  ...
  A = Jacobian(f(...)) # A comes from some Jacobian
  K = linear_solver(A) # Build the factorization object
  ... # Some calculations
  @into! x1 = K\b1 # First linear solve, default settings from factorization object
  ... # use x1 a bit
  solve!(x,K,b2,MaxIterations(10)) # Suppose the second linear solve just needs an approximation, stop after 10 iterations
  ... # use x2 a bit
  solve!(x,K,b3,Transpose(true)) # third linear solve uses the transpose of A
  ... # Do some stuff to end
end

```

The same could probably be achieved by modifying `A`, but this would be trickier if the solutions could happen in parallel. In general, I think it’s prudent to have a provision for passing parameters related to the application of `A` for a specific RHS but that are not directly a property of `A` itself.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 12, 2016, 7:03pm UTC](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699/27 "2016-12-12T19:03:23Z")

</div>

> [@barche](#):
>
> You could set e.g. tolerance at each step:

What about just requiring that every factorization object implements an `update!` function?

> [@ChrisRackauckas](#):
>
> function OneStep(f,…,linear\_solver=lufact,update\_args\_1=nothing,update\_args\_2=nothing)  
> …  
> A = Jacobian(f(…)) # A comes from some Jacobian  
> K = linear\_solver(A) # Build the factorization object  
> … # Some calculations  
> @into! x1 = K\b1 # First linear solve  
> … # use x1 a bit  
> update!(K,update\_args\_1…)  
> @into! x2 = K\b2 # Second linear solve  
> … # use x2 a bit  
> update!(K,update\_args\_2…)  
> @into! x3 = K\b3 # third linear solve  
> … # Do some stuff to end  
> end

The only thing required to make this work well would be a PR to Base convincing them that a blank dispatch for `update!` would make this code much easier to write (and maybe make a new spot on the “Interfaces” page of the documentation about this interface).

---

<div class="post-metadata">

**Author:** ![barche](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/barche/32/79_2.png) [@barche](https://discourse.julialang.org/u/barche)\
**Post date:** [December 12, 2016, 7:44pm UTC](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699/28 "2016-12-12T19:44:40Z")

</div>

(I Realize I wrote `A` instead of `K` in my previous post, will correct after sending this)

This does not address the case (a bit far fetched, admittedly, and not something I require in my code) where the solves would happen in parallel: different parameterizations of `K` would require making copies. Also, conceptually things like the number of iterations and the tolerance aren’t really a property of the factorization object, but of a specific application of that object, so I feel it makes sense to be able to set these parameters when they matter (allowing of course for defaults to be set in `K`, to allow the `\` variant to work).

Regarding `update!`: How about a `parameters` function that gets parameters that can be manipulated using `setindex!`:

```julia
params = parameters(K)
setindex!(params, 1e-6, :tolerance)

```

The name `update!` seems a bit ambiguous, it could also be interpreted as “I changed the matrix, update the factorization”, for example.

Also, is it necessary to have all the interface for this in Base? Couldn’t a package describing the interface be created? It should then be easier to make changes if necessary.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 12, 2016, 7:55pm UTC](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699/29 "2016-12-12T19:55:24Z")

</div>

> [@barche](#):
>
> Also, is it necessary to have all the interface for this in Base? Couldn’t a package describing the interface be created? It should then be easier to make changes if necessary.

If we want to use \ and `A_ldiv_B!` (which I still think is sensible), we would only need a minimal amount of Base compatibility (i.e. just ways for it to ignore the parameters instead of throwing an error when they’re there).

---

<div class="post-metadata">

**Author:** ![barche](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/barche/32/79_2.png) [@barche](https://discourse.julialang.org/u/barche)\
**Post date:** [December 12, 2016, 8:07pm UTC](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699/30 "2016-12-12T20:07:57Z")

</div>

> [@ChrisRackauckas](#):
>
> If we want to use \ and A\_ldiv\_B! (which I still think is sensible), we would only need a minimal amount of Base compatibility (i.e. just ways for it to ignore the parameters instead of throwing an error when they’re there).

Yes, the parameters should always be optional, with defaults configurable in `K`. I agree that \ and `A_ldiv_B!` should stay, of course, but the `solve!` interface would be a nice addition to add a little more flexibility. I must admit I never would have thought to look for a function named `A_ldiv_B!`, would it be a totally crazy idea to rename that to `solve!`?

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [December 12, 2016, 9:32pm UTC](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699/31 "2016-12-12T21:32:27Z")

</div>

> [@barche](#):
>
> A\_ldiv\_B!, would it be a totally crazy idea to rename that to solve!?

I’m not a huge fan of the naming scheme, but it’s consistent with `A_mul_B!`, `Ac_mul_B!`, etc. The notation `solve!` on its own is not as versatile.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 13, 2016, 1:50am UTC](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699/32 "2016-12-13T01:50:09Z")

</div>

In parallel to the discussion on parameters, I think we should bring up the discussion of linear maps. Not every linear system that is to be solved is a matrix, in many cases one might want to solve with a linear operator or a linear map. Generalizing it like this computationally good because things like matrix-free methods can be much faster and/or take much less memory.

IterativeSolvers.jl seems to want to go this direction by using LinearMaps.jl

> **[GitHub - JuliaLinearAlgebra/LinearMaps.jl: A Julia package for defining and...](https://github.com/JuliaLinearAlgebra/LinearMaps.jl)**
>
> A Julia package for defining and working with linear maps, also known as linear transformations or linear operators acting on vectors. The only requirement for a LinearMap is that it can act on a v...

I think it would be nice if more packages try to support at least something like linear maps. If I understand correctly, a `WrappedMap` is just a matrix wrapped as a linear map, so I assume that the package handles the cases to make this work anywhere a matrix \ would work (does it?).

Then the push would really be about trying to get other libraries to support something like the `FunctionMap` where you define a function for `A*x`. It’s clear that iterative solvers “can” work with this, but what packages can be made compatible with something like this? Do things like CUSOLVER and PETSc strictly require a matrix?

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [December 13, 2016, 4:52am UTC](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699/33 "2016-12-13T04:52:44Z")

</div>

I’m confused why LinearMaps.jl doesn’t do the following:

```julia
abstract AbstractLinearMap{T} <: AbstractMatrix{T}

getindex(A::AbstractLinearMap,k::Integer,j::Integer) = (A*[zeros(j-1);1;zeros(size(A,2)-j)])[k]

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 23, 2017, 7:34pm UTC](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699/34 "2017-01-23T19:34:54Z")

</div>

An update on this: the factorization-type style has been implemented in DifferentialEquations.jl. I tested it out, it works really nicely as a cross package solution. While it’s only currently enabled on master, docs are already up which show how to specify it:

[http://docs.juliadiffeq.org/latest/features/linear\_nonlinear.html](http://docs.juliadiffeq.org/latest/features/linear_nonlinear.html)

I’ll see if I get around to putting some PRs around to make this work with more packages. All that needs to happen is the type interface like PETSc.jl needs to be created. There’s a thread going in IterativeSolvers.jl for this:

[https://github.com/JuliaMath/IterativeSolvers.jl/issues/99](https://github.com/JuliaMath/IterativeSolvers.jl/issues/99)

I’ll also be poking around for packages which should use this interface, like NLsolve.jl.

[https://github.com/EconForge/NLsolve.jl/issues/83](https://github.com/EconForge/NLsolve.jl/issues/83)

In total, it’s really painless to setup but is a really cool feature for Julia. Thanks everyone for the ideas!

---

<div class="post-metadata">

**Author:** ![JaredCrean2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jaredcrean2/32/3574_2.png) [@JaredCrean2](https://discourse.julialang.org/u/JaredCrean2)\
**Post date:** [January 25, 2017, 3:04am UTC](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699/35 "2017-01-25T03:04:50Z")

</div>

@ChrisRackauckas Re. your earlier question, PETSc requires a matrix object. Matrix objects can have different formats, including various flavors of sparse storage formats and MatShell, which doesn’t require an explicit matrix but instead requires a callback function that performs matrix-vector products.

[Previous page](https://discourse.julialang.org/t/unified-interface-for-linear-solving/699.md?page=1)