A_ldiv_B!(a, b) writes over b or not depending on type of a

[KristofferC]

Is this intended?

julia> a = rand(3,3); b = rand(3)
3-element Array{Float64,1}:
 0.95134 
 0.43887 
 0.719551

julia> A_ldiv_B!(sparse(a), b) 
3-element Array{Float64,1}:
 0.455667
 0.712479
 0.326213

julia> b # Note: b is not overwritten
3-element Array{Float64,1}:
 0.95134 
 0.43887 
 0.719551

julia> A_ldiv_B!(lufact(a), b)
3-element Array{Float64,1}:
 0.455667
 0.712479
 0.326213

julia> b # For this, b is overwritten
3-element Array{Float64,1}:
 0.455667
 0.712479
 0.326213

[dmbates]

It's a misnomer in the Umfpack interface, which is what is used for solving square, non-Hermitian, sparse systems. The A_ldiv_B! method calls the internal solve function in base/linalg/umfpack.jl which does not overwrite the right hand side but instead returns a freshly allocated vector.

Whether or not this is a bug depends on how strongly you interpret the ! in the name A_ldiv_B!, etc. Does it mean that it can overwrite the argument B or does it mean that it must overwrite the argument B?

[dmbates]

I should add that, as far as I know, there is no mutating solver in Umfpack, so the only way you could get mutating behavior in Julia is to allocate a vector or matrix for the solution, solve the system, then copy the solution back into the right hand side.

[tkelman]

I think a no-method error would be preferable in cases where a mutating operation is not available, rather than trying to mimic an in-place operation by allocating and copying.

[toivoh]

I'm not sure. Let's say you wrote an algorithm using the updating versions for increased efficiency. Wouldn't you still want it to work for other cases, even though you can't get at the same efficency? Btw if there is no mutating solve in umfpack, maybe Tim Davis would be sympathetic to adding one? But I guess it's not as useful if it's only in new releases.

[tkelman]

Having the in-place version perform slower and allocate [as much or] more memory than the standard case defeats the purpose of having it at all. If you want code to be generic and not every type can implement an updating operation, then you should probably be using the conventional operation. If you want the highest possible performance, your code would need to be aware of the types it's operating on, use in-place for types where it is faster, or the standard form when an in-place operation is not available.

Yes it would be ideal if every single type supported in-place operations and they would always be the higher-performance option, then the copying operation would be always be a trivial generic wrapper around the in-place version. But this is not currently the case.

edit: The current situation of having an in-place operator that isn't actually in-place is semantically confusing, but at least doesn't incur a performance cost due to extra copying. So maybe it's acceptable for now?

[nalimilan]

@tkelman That's a tough call. You may want your algorithm to be as fast as possible for types that support in-place operations, and still work (even if less efficiently) for other types. The alternatives are 1) the function fails plainly when given a type that doesn't support in-place operation, or 2) the developer checks for the existence of the in-place method for the type, and falls back to the copying more generic function.

But I agree that A_*B! functions imply in their name that the operation is in-place, and that silently making a copy with some types can be confusing. Maybe with an improved syntax for in-place operations (#249), falling back to an allocating version when needed would be more acceptable.

[toivoh]

The semantics would of course still have to be to do the update in place.

[KristofferC]

Personally, I expect some homogeneity in the A_op_B methods. It feels dangerous if you try to write a general method for matrices and you get different outputs depending on the representation you chose for the matrix.

I think that there are two possible choices:

• Don't have a mutating method for the types where there is no in-place solver.
• Manually do the mutation of, in this case, b.

Keeping the current situation feels like it can lead to very subtle bugs in user code.

[tkelman]

Manual mutation and incurring the associated performance penalty might be preferable to a no-method error, but then you have a more subtle performance discrepancy that depends on the input types.

[toivoh]

We do provide generic fallbacks for some AbstractArray methods. I think this could be seen to be in a similar spirit.

Btw had anyone timed this to see how much efficiency is actually lost in this case?

[hayd]

I always read ! as "may modify some/all arguments" rather than "inplace", specifically in this case I've no reason to think that a remains fixed unless I read the documentation (or source). For code clarity I'd always write:

b = A_ldiv_B!(a, b)  # or perhaps I'd assign it to a new variable

Similarly you don't really lose anything by writing:

s = sort!(s)

but the documentation says clearly you don't need to.

[toivoh]

In general I agree that the ! doesn't mean that a function has to modify to its arguments, but for in place versions of operations like A_mul_B!, I think it has come to be expected. It seems like an invitation to subtle bugs if we make the mutation optional.

I also seem to remember from a previous discussion that always assigning the return value back to the corresponding variable might make some optimizations harder for the compiler.

[KristofferC]

Not sure how much it matters, but in Python it seems that methods that operate inplace usually do not return anything, for example list.sort().

If this is good or not, I am not sure. Many times students write something similar to:

def sort_thingy(thingy)
    return thingy.sort()

and get confused when they get None popping up here and there.

[timholy]

Agreed that A_ldiv_B! needs to be consistent: might mutate is sort of like might copy in reshape, and it seems we all agree that needs to change (#10507).

@KristofferC, I think it's better to return something, that way you can do mean(sqrt!(x)).

[KristofferC]

@timholy I agree that is is more comfortable to return the argument since you can chain expressions like you mentioned. @toivoh mentioned that it might come at a performance cost. Anyone know more about that?

[pao]

mentioned that it might come at a performance cost

For the explicit reassignment version only, it sounds like, but...

I'm a bit worried about mean(sqrt!(x)), personally--it's probably slower than the equivalent loop, and mutates x as a side effect but not obviously so because the sqrt!() is buried in the expression. Consenting adults and all that, but I wouldn't encourage that style. (After all, the standard for research on human subjects is informed consent. And programmers are still humans.)

[toivoh]

The performance cost I was talking about was for assigning the return value of a mutating function back to the original variable - the compiler might not always be able to realize that it is a noop. But I don't think we should encourage that style anyway. Python is kind of purist when it comes to returning None (python nothing) from mutating functions, to avoid hiding the side effects. But in practice, chaining can be pretty useful. I think our case is slightly easier at least thanks to the visually distinguishable exclamation mark.

[matthieugomez]

+1. It's a source of bug. It also makes it harder to understand which operations really allocate in place.

[KristofferC]

I close this now, see #13496