inv() of one-element Array

[pao]

Reported by Eduardo H on julia-users https://groups.google.com/d/msg/julia-users/YIY9oXolMMU/0gdpaNnChsoJ

julia> x = [1:3]
3-element Int64 Array:
 1
 2
 3

julia> inv(x'x)
ERROR: no method inv(Array{Int64,1},)

Should probably be [0.0714286].

[andreasnoack]

I tried to browse the mailing list but couldn't find the argument for v'v returning a Vector. The calculation dispatches to Ac_mul_B and not transpose(A)*B so we can choose the return method freely. I would prefer v'v to be a scalar but have no problem with it being a 1x1 Matrix. However, I find it a bit weird that it is a Vector.

If v'v continues to be a Vector, should the inverse then be a scalar, vector or matrix?

[kmsquire]

+1 on returning a scalar for v'v

[kmsquire]

Or any vector-vector inner product!

[StefanKarpinski]

The correct behavior depends whether we view v'*w as a single '* "inner product" operation or as just a shortcut for computing the product after taking the transpose of the left operand. We have been thinking of it as simply a shortcut, in which case the transpose creates a matrix and the product is a matrix-vector product, which should produce a single-element vector. However, maybe things like '*should be considered first-class operations.

Another "way out" would be to make transpose lazy, so that v' can be a value of type Transpose{Vector}. This has a lot of benefits, including allowing us to remove the special parsing of v'*w and such altogether. Another benefit is that v'' can produce a vector instead of producing a column matrix.

[toivoh]

+1 for Transpose{Vector}! Should work even better with Transpose as an immutable type.

[kmsquire]

@StefanKarpinski, that makes sense--thanks for the explanation. I hadn't thought through the implications (need my coffee...) I also like the lazy transpose idea. It's been a while since I touched LAPACK directly, but from other discussions here and from what little I remember, LAPACK should actually support Transposed matrices directly, right? Assuming someone took the time to write the code in Julia.

[StefanKarpinski]

As @toivoh pointed out, you need an immutable Vector type (possibly immutable by convention, wrapping a mutable array like our strings do), for this to work well. Otherwise you still end up needing to make a copy when you do the "lazy" transpose. Which isn't really very lazy.

[andreasnoack]

@kmsquire Actually some functionality for transposed matrices in LAPACK and BLAS is already in use through the Ac_mul_B type functions which are called when you write A'B.

I like the Transpose type and can volunteer to go through the linalg code if we get it.

[lindahua]

+1 for the Transpose type. Perhaps, we may also have ConjugateTranspose for complex vectors/matrices.

[StefanKarpinski]

The Transpose type isn't really a win unless we adopt a distinction between mutable arrays and immutable tensors wrapping mutable arrays, making them immutable by convention. Otherwise we would end up doing more copying than we do now.

[toivoh]

We would need a way to force the transpose when you want an actual Array
as a result. Is that what you mean?

Beside that, a Transpose object wouldn't have to act like an immutable
array, but could be a transposed view on the underlying data. E.g.
assigning into it would permute the indices and then assign on the
underlying array. But I'm not sure if that would be a gain or an efficiency
loss.

When I was talking about immutable types, it was only as an implementation
detail to mark the outstanding transpose operation. @StefanKarpinski: Do
you think that the Transpose object should provide a const view to the
underlying array?

[toivoh]

Simple ways to force the transpose would probably be convert(Array, x') (which might be called implicitly in a number of places?), and [x'].

[JeffBezanson]

Of course a 1x3 matrix times a 3-vector gives a 1-vector. But IIUC the transpose of a vector is not necessarily a 1xN matrix. So it seems defensible for (transposed vector)*vector to be different than (1xN matrix)*vector (and give a scalar result).

It would also be ok for inv of a 1-vector to work I think.

[toivoh]

I think that it might make sense to have inv work on 1-vectors, but then it should do so by extracting the single scalar out of the vector, and calling inv on that. The reason that it would make sense in the first place is that using inv asserts that the one dimension that the vector has was accidental. (Ok, it could alternatively convert the vector to a 0d array)

[farshid-vahid]

x=[1:2]
inv(x'x)

still returns an error. Why is this issue closed?

[KristofferC]

0.6-rc2:

julia> x=[1:2]
1-element Array{UnitRange{Int64},1}:
 1:2

julia> inv(x'x)
0.2

[farshid-vahid]

Thanks. I was using 0.5.2. Cheers.