Support `I * X = X` for Arrays, i.e., UniformScaling as a true identity operator

[JeffFessler]

The docstring for UniformScaling describes it as an "identity operator"
yet the support multiplication works only for a Vector or Matrix,
not for a general Array, as illustrated by this MWE:

using LinearAlgebra: I
5I * ones(2,3) # works of course
5I * ones(2,3,4) # fails with MethodError

The culprit is this one line of source code:

*(J::UniformScaling, A::AbstractVecOrMat) = J.λ*A

LinearAlgebra.jl/src/uniformscaling.jl

Line 279 in 27b0f13

*(J::UniformScaling, A::AbstractVecOrMat) = J.λ*A

that limits A to be a matrix or vector.

My feature request is to generalize this to

*(J::UniformScaling, A::AbstractArray) = J.λ*A

With this tiny change, I really become an "identity operator" rather than merely an identity "matrix."

For context, I've put this more general method in one of my packages for years, but Aqua.jl tells me that it is type piracy, so it would be better to have this generality here in LinearAlgebra.
I will make a PR.

[stevengj]

On the one hand, multi-dimensional arrays do form a vector space in Julia, so it could make sense to define linear operators on them, and of course scaling is one such operation.

On the other hand, we don't define any generic linear operations on multidimensional arrays, and it's not clear how useful this would be in isolation? If you just want to scale an array you can use the vector-space operation of multiplying it by a scalar.

Do you have a concrete application in mind? It is hard to discuss a feature in the abstract.

What do packages that define generic linear operations or multi-linear algebra do now for the identity?

[JeffFessler]

Thanks for the quick feedback.
To try to answer your question about other packages, here are two (closely related) examples:

https://github.com/JuliaLinearAlgebra/LinearMaps.jl/blob/d34be208afac772014d522155eb1722c5f9a4143/src/uniformscalingmap.jl#L33

https://github.com/JeffFessler/LinearMapsAA.jl/blob/97b2518139753a238bec1bb84b6d3a7100c1a35f/src/multiply.jl#L21

In these operator packages it is helpful to have support for I because more complicated operators can be constructed from I, e.g. (I + A) * B etc., and then when multiplying that composite operator by a tensor eventually there will be an I * X step. It is correct to write I + A to form a (lazy) operator, whereas constant + A is incorrect (as you know of course).

we don't define any generic linear operations on multidimensional array

That's true, but we do support dot(A,B) for multidimensional arrays.
BTW, for consistency, if we do this we should also support X * I = X.

[stevengj]

BTW, for consistency, if we do this we should also support X * I = X.

It's not clear to me what right-multiplication means here?

[JeffFessler]

BTW, for consistency, if we do this we should also support X * I = X.

It's not clear to me what right-multiplication means here?

Oops, never mind about right multiply 😄 It makes sense for matrix algebra but not for general linear algebra. Sorry!

[dkarrasch]

I'm slightly confused. In LinearMaps.jl, we have the UniformScalingMap type, so in your composite example it should happen that at some point a UniformScalingMap hits a tensor, which you could define as scalar multiplication, without relying on LinearAlgebra.UniformScaling? If you remove your type piracy definition, what's the stack trace that comes up?

[JeffFessler]

There's no problem with LinearMaps AFAIK. I just brought that up to answer @stevengj's question by giving examples of places where identity-type operators are used in practice.

@dkarrasch, leaving aside type piracy and any specific packages, what is your view about extending LinearAlgebra.UniformScaling so that it truly is an identity operator as advertised in its docstring, instead of just an identity matrix? It is a tiny change to the source code that is 100% backwards compatible with existing behavior. It adds linear algebra functionality without any performance penalty (unless I'm missing something, which is certainly possible).