check sizes of arguments in dot; fixes #28617

[ranocha]

As discussed in #28617 and https://discourse.julialang.org/t/efficient-trace-of-product-of-matrices/13313/12, dot(A::AbstractArray, B::AbstractArray) should check whether size(A) == size(B) instead of length(A) == length(B).

[andreasnoack]

I think the check should actually be axes(A) == axes(B) to handle non-standard indices correctly.

[ranocha]

I think the check should actually be axes(A) == axes(B) to handle non-standard indices correctly.

I thought custom indices should only be some kind of simplification for special purposes. Thus, I would not forbid having two arrays of the same size but with different axes and calculating the dot product of them. Could you please explain why we should disallow something as

julia> using LinearAlgebra, OffsetArrays
                                                                                    
julia> A = [1 2; 3 4]; B = OffsetArray(A, (-1,-1))
OffsetArray(::Array{Int64,2}, 0:1, 0:1) with eltype Int64 with indices 0:1×0:1:
 1  2
 3  4

julia> dot(A, A)
30

julia> dot(A, B)
30

[andreasnoack]

The idea is that the indices of A and B should match. So for Array, we'd disallow dot(A, B) for A = randn(4) and B = randn(2,2) because vec(CartesianIndices(A)) != vec(CartesianIndices(B)). The generalization of this requirement to AbstractArray is to require axes(A) == axes(B). I'd also expect this to hold when working with OffsetArrays but I might, of course, be wrong. Eventually, I'd also think it would be reasonable to require that axes(A, 2) == axes(B, 1) in the matrix multiplication A*B.

[ranocha]

Thank you for the explanation. If this is the basic idea behind the checks (and A*B will work if axes(A, 2) == axes(B, 1)), there should be checks using axes, of course. I will update the PR accordingly.

[garrison]

"of the same size" -> "with the same axes" ? The "equal sizes" language is also used a few lines below...

[ranocha]

I will change that.

[garrison]

Is axes purposely not used here for some reason?

[ranocha]

I thought SparseMatrixCSC uses definitely classical indexing as in line 208. Thus, I've just used size instead of axes. Should I change that?

[jebej]

For now, the SparseArrays code has not been modified to support arbitrary indices.
I suppose it would be generally less useful than for dense arrays.

[ranocha]

We should note somewhere that these changes are breaking. Since the NEWS.md file seems to be reserved for the official release of Julia v1.0, I don't know where I should add a comment. Could you help me, please?

[stevengj]

I guess the question is whether we consider this to be a bugfix or a change in the documented API.

[ranocha]

Following the discussion above, it seems that most people see this as a bugfix. Should I modify this PR in some way before it can be merged?

[KristofferC]

This is not a simple bug fix.

For any iterable containers x and y (including arrays of any dimension) of numbers (or any element
type for which dot is defined), compute the dot product (or inner product or scalar product), i.e. the
sum of dot(x[i],y[i]), as if they were vectors.

It works exactly like documented. We can't label things bugfixes just because we change what we like about the behavior. Also, master is currently closed for breaking changes.

[ranocha]

Okay. Do I have to anything now in that case?

[KristofferC]

This LGTM so I think this can rest for a while until we open up master or implement versioning in for the stdlibs.

[garrison]

Bump. Any chance this qualifies as a "minor change"?

[StefanKarpinski]

I've marked it as a "minor change" and put the "triage" label on it and 1.1 milestone; we can discuss whether to include it on Thursday's triage call.

[JeffBezanson]

I'm not sure this is useful at all --- it just adds errors to more cases. Very much not sure it's worth the potential breakage.

[JeffBezanson]

Also not obvious to me whether comparing axes or just sizes is the right thing.

[StefanKarpinski]

The general triage sentiment is that by default we need to lean heavily in favor of not breaking code. So this could happen but the burden is on those proposing this that this is (a) the right thing to do and (b) that it won't break any packages or applications. Running PkgEval on this branch would be a good first step to inform that decision—if nothing breaks, we can consider it.

[JeffBezanson]

Triage says 👍 but we need a PkgEval run.

[marius311]

I'm for the behavior in this PR. Also, It should be marked as fixing #32395 (which I opened).

[andreasnoack]

@ararslan Would you be able to run PkgEval here?

[ararslan]

I won't likely be able to get to it for a while.

[maleadt]

@nanosoldier runtests(ALL, vs = ":master")

[nanosoldier]

Your test job has completed - possible new issues were detected. A full report can be found here. cc @maleadt

[JeffBezanson]

Ok, this does indeed cause some package test failures. That matches my intuition that, given that we allow passing arbitrary iterators to dot, being permissive here can be useful.

[andreasnoack]

I just read through the error logs here and I think we have drawn the wrong conclusion from them. Most of the failures are unrelated to this PR. They are either segfaults or look intermittent. Then three of them fail because of an unfortunate sum method in StatsBase that was removed in JuliaStats/StatsBase.jl#526. Finally, it looks like there are just a few cases where somebody ends up dotting a vector and a one-column matrix. The latter might be something that people actually find convenient but let's rebase one and rerun PkgEval to get an updated view on this.

[maleadt]

Third time's the charm (and comes with a prevind fix):

@nanosoldier runtests(ALL, vs = ":master")

[nanosoldier]

Your package evaluation job has completed - possible new issues were detected. A full report can be found here.

[andreasnoack]

I've tried to summarize the reason for the failures in the table below. No packages rely on dotting arrays of different shape except for singleton dimensions. A few packages rely on dotting one column matrix with a vector and a few packages dots a one row matrix with a vector. There are also some numerical differences which might be because some calls now use our dot instead of the dot from BLAS. We can revert that but I think the tests might just be too strict and that it would be better to go through the packages and adjust the tests.

I think it's unfortunate that we don't provide an error in the motivating example in #28617 but maybe it's not a sufficiently common case to warrant the the effort required for introducing the check.

Package	Reason	src or test
AeroMDAO	One column matrix	src
CompressiveLearning	One column matrix	src
CopEnt	Unrelated	N/A
EnhancedGJK	One column matrix	src
ExponentialUtilities	Unrelated	N/A
FaultDetectionTools	Unrelated	N/A
FilesystemDatastructures	Unrelated	N/A
GridArrays	Unrelated	N/A
ImageDistances	Unrelated	N/A
ImmersedLayers	Marginally different results	N/A
IndependentComponentAnalysis	Unrelated	N/A
KissMCMC	Marginally different results	N/A
LearningHorse	One row(?!?) matrix	src
MultinomialRegression	Marginally different results	N/A
MutableArithmetics	One column matrix	test
OpSel	One column matrix	src
QuantumTomography	Unrelated	N/A
RegressionDiscontinuity	One row(?!?) matrix	src
SolverTools	Marginally different results	N/A
Symbolic	Unrelated	N/A
TensorCore	One row(?!?) matrix	Test
UnbalancedOptimalTransport	One column matrix	src
WaveFD	Marginally different results	N/A
Widgets	Unrelated	N/A