Improve matrix multiplication and solvers

[bardo84]

See

https://discourse.julialang.org/t/julia-vs-zig-surprise/74540/11
The Linear Algebra Mapping Problem. Current State of Linear Algebra Languages and Libraries

[antoine-levitt]

Would be good to be clear what exactly is missing. Julia can definitely do diagonal/triangular matrix multiplies, it just wants you to use the appropriate types (ie Diagonal instead of plain matrix). It doesn't try to read off properties from the numerical values of the matrix (eg it won't do anything fancy with diagm), which I think is completely reasonable.

[pauldj]

Hi Antoine, I am one of the authors of the paper. While it is perfectly reasonable to expect users "to use the appropriate types" (as you say), there are many situations when this is not possible. For instance, when the matrices are the result of intermediate computations.

As a simple example, given an input matrix X, the result of X * X^T is symmetric. If a user writes (X * X^T) * B, it would be beneficial to take advantage of symmetry. Notice that this is only the simplest of the examples, and could potentially be handled by having the kernel that computes X * X^T (SYRK) return also a flag ("symmetric") for the output matrix. In general this is not possible. One extremely common subexpression that arises in many practical cases is (X * D * X^T). In this case, returning a flag is not a viable solution.

In general, there are tons of situations in which one operates with matrices without knowing their properties. Some of our experiments test the behaviour of the different languages in these conditions.

[hbarthels]

Hi Antoine,
I'm another one of the authors. I'm not sure if we cover this in the paper, but in practice it can also happen that matrices have more than one property. As an example, a matrix could be diagonal and symmetric positive definite (SPD) at the same time. In this case, this just means that all the diagonal elements are positive, but of course there are other combinations.

[mcabbott]

About X * X^T, I believe the reason these don't return a Symmetric is that the type system doesn't know that X === X i.e. that the same matrix appears twice. But the syntax does know this, so one easy step would be to write a macro which does this to whole blocks of code. A bigger step would be to build that into the parser, like powers:

julia> 2^-1  # magic literal_pow
0.5

julia> 2^(0-1)  # simple ^ can't do that, as it would be type-unstable
ERROR: DomainError with -1:
Cannot raise an integer x to a negative power -1.
Make x or -1 a float by adding a zero decimal (e.g., 2.0^-1 or 2^-1.0 instead of 2^-1), or write 1/x^1, float(x)^-1, x^float(-1) or (x//1)^-1

The same goes for X * D * X', which is parsed as one operation, before being broken up.

Agree these wrappers don't compose easily. It doesn't sound crazy to encode positive-definite somehow, I wonder if this was discussed? Closest I can find is https://github.com/JuliaLang/julia/issues/20645 where the example is:

julia> m, n = 2, 3;

julia> A = randn(m,n) + im*randn(m,n);

julia> M = A * Diagonal(rand(n)) * A'
2×2 Matrix{ComplexF64}:
   5.57309-2.22045e-16im  -0.974748-1.35619im
 -0.974748+1.35619im       0.948511+0.0im

julia> eigvals(M)
2-element Vector{ComplexF64}:
    2.2811235643286 - 3.3197837198443943e-16im
 15.371828833164777 + 3.319783719844396e-16im

julia> Hermitian(M)  # hides the imaginary part of diagonal entries
2×2 Hermitian{ComplexF64, Matrix{ComplexF64}}:
 10.5119+0.0im      3.72908+5.10833im
 3.72908-5.10833im  7.14107+0.0im

julia> eigvals(Hermitian(M))
2-element Vector{Float64}:
  2.281123564328599
 15.371828833164772

[Roger-luo]

Just a side comment about the symmetry property I think this type of property be inferred by abstract interpreter or some intermediate IR for arrays (like MLIR) in principle? So that one doesn't need to write it manually. E.g the above example can be matched with A * A' pattern. We do this type of thing more aggressively for quantum circuits since it happens almost everywhere. So I think the more general case for matrices can also be handled by a proper compiler transform.

[oscardssmith]

I believe Julia is currently missing optimized Symmetric*Symmetric.

[bardo84]

I have contacted the authors of the paper for input.

[mcabbott]

Table 7: Experiment #6 in the paper is also wrong, for Julia 1.7.

Edit: at last I found it, they specify Julia 1.5.2, thus too old for this.

Their article is "[v2] Sun, 5 Sep 2021" without comment on what changed, perhaps they did not re-run experiments.

[pauldj]

This is very nice! Maybe our paper had an impact :-)

[ViralBShah]

May I suggest that we use discourse to zero in on the problems and then open specific issues for each case that are actionable? If folks agree, we can close this issue.

[DilumAluthge]

Alternatively, if you like, we can convert this issue to a GitHub Discussion - just click on the "convert to discussion link", which will close the issue and convert it to a discussion.

[mkitti]

I believe my colleague, @bjarthur, is working on symmetric matrix multiplication. We currently have a working private repository for this. I will ping him about this discussion.

[bjarthur]

in the OP's linked discourse thread i have just described my PackedArrays.jl [EDIT: SymmetricFormats.jl] package.

@DilumAluthge i don't think another discourse thread is necessary. and thanks @mkitti for the heads up.