Add dot method for symmetric matrices

[goggle]

This PR adds a dot method for symmetric matrices, which makes the calculation of the dot product of two symmetric matrices significantly faster (see benchmarks below). It intends to fix #32730.

To perform the benchmarks, this code was used:

using LinearAlgebra
using BenchmarkTools

n = 10_000
X = randn(n, n)
Y = randn(n, n)

A = Symmetric(X)  # or A = Symmetric(X, :L)
B = Symmetric(Y)  # or B = Symmetric(Y, :L)
@btime dot(A, B)

These are the results:

A.uplo	B.uplo	Julia 1.1	PR
'U'	'U'	1.4 s	72 ms
'U'	'L'	1.33 s	537 ms
'L'	'U'	1.33 s	538 ms
'L'	'L'	1.4 s	72 ms

[dkarrasch]

This looks good, but I guess a few changes are required for block matrices, taking into account that dot acts recursively.

[goggle]

I have done the requested changes by @dkarrasch.

I don't quite understand the possible problem with block matrices. I didn't find them in the Julia stdlib, so I guess you mean the package BlockArrays.jl @dkarrasch? If you think there is an issue, feel free to elaborate.

[dkarrasch]

No, I meant matrices of matrices. Like

X = [rand(ComplexF64, 10, 10) for _ in 1:3, _ in 1:3]
dot(Symmetric(X), Symmetric(X)) # should give a real, non-negative number

[antoine-levitt]

Could this also support Hermitian?

I think it'd be good to add more comprehensive tests, including the four code paths (L/R), complex and block matrices: mistakes in functions like these might not be trivial to spot, and could have very bad consequences.

[goggle]

@dkarrasch Thanks for the clarification, I'm going to check this!

@antoine-levitt Sure, I can add a similar method for Hermitian, but I would prefer to do this in a followup-PR. I agree on improving the tests.

[chriscoey]

@goggle thanks for this! It speeds up my code.
I agree with@antoine-levitt's points. Hermitian should be an easy generalization.

[goggle]

There is indeed an issue with block matrices.
Let's consider consider this example:

using LinearAlgebra
a = [1 0 ; 1 0]
b = [-1 1; 0 1]
c = [-1 -1 ; 0 0]
d = [-1 2 ; -1 0]
A = [[a] [c] ; [b] [d]]
symA = Symmetric(A)

symA is now a matrix consisting of four (2x2)-blocks, such that symA is symmetric when considered as a 4x4 matrix. In the symA.data fields we have the origin data (the arrays a, b, c and d), so my method produces a wrong answer: dot(symA, symA) = 14 on Julia 1.1.1 vs. dot(symA, symA) = 12 on 738b89a.

I will replace all the A.data and B.data by A and B. This should produce correct results, but might be slightly slower, so I will also rerun the benchmark.

[goggle]

Alright, so I avoid to access the data directly, which should solve the issues with symmetric block matrices. As expected, it makes it slightly slower than the previous version. Here are the updated benchmarks:

A.uplo	B.uplo	Julia 1.1	bae0753	dcc1aae
'U'	'U'	1.4 s	72 ms	109 ms
'U'	'L'	1.33 s	537 ms	602 ms
'L'	'U'	1.33 s	538 ms	608 ms
'L'	'L'	1.4 s	72 ms	100 ms

I will extend the tests tomorrow.

[goggle]

It took me a while to resolve many different issues.
I have restricted this specialized dot method to symmetric matrices that have a Number type in their entries. The reason for this is, that the suggested algorithm does not work for symmetric block matrices: The blocks itself are not symmetric, only the matrix as a whole.
There are also more tests now.
It also uses the data field directly again, so the benchmarks from #32827 (comment) are valid again.

This PR can now be reviewed again (@dkarrasch).

[dkarrasch]

I don’t see any need to restrict to numbers. The remaining issues for nested arrays are easily fixed. What was missing is using the implicit symmetry of diagonal elements, which may not be already contained in the data field. For the Hermitian version of this (not sure if we would also add mixed versions), you will need to replace the off-diagonal contributions by

dotprod += 2real(dot(A.data[i, j], B.data[i, j]))

otherwise this should be the same.

[goggle]

Thanks for these detailed comments @dkarrasch, this helps a lot. I would also prefer to have a generic implementation dot(A::Symmetric, B::Symmetric) for the dot product of symmetric matrices, but I have run into some issues that still persist.

Here is the modified dot method with your suggestions from #32827 (comment) and #32827 (comment) applied:

using LinearAlgebra

function mdot(A::Symmetric, B::Symmetric)
    n = size(A, 2)
    if n != size(B, 2)
        throw(DimensionMismatch("A has dimensions $(size(A)) but B has dimensions $(size(B))"))
    end

    dotprod = zero(dot(first(A), first(B)))
    @inbounds if A.uplo == 'U' && B.uplo == 'U'
        for j in 1:n
            for i in 1:(j - 1)
                dotprod += 2 * dot(A.data[i, j], B.data[i, j])
            end
            dotprod += dot(A[j, j], B[j, j])
        end
    elseif A.uplo == 'L' && B.uplo == 'L'
        for j in 1:n
            dotprod += dot(A[j, j], B[j, j])
            for i in (j + 1):n
                dotprod += 2 * dot(A.data[i, j], B.data[i, j])
            end
        end
    elseif A.uplo == 'U' && B.uplo == 'L'
        for j in 1:n
            for i in 1:(j - 1)
                dotprod += 2 * dot(A.data[i, j], B.data[j, i])
            end
            dotprod += dot(A[j, j], B[j, j])
        end
    elseif A.uplo == 'L' && B.uplo == 'U'
        for j in 1:n
            dotprod += dot(A[j, j], B[j, j])
            for i in (j + 1):n
                dotprod += 2 * dot(A.data[i, j], B.data[j, i])
            end
        end
    end
    return dotprod
end

I call it mdot here. Let's take the block matrix from #32827 (comment):

a = [1 0 ; 1 0]
b = [-1 1; 0 1]
c = [-1 -1 ; 0 0]
d = [-1 2 ; -1 0]
A = [[a] [c] ; [b] [d]]
sau = Symmetric(A, :U)
sal = Symmetric(A, :L)

Now I get a wrong result with mdot:

julia> mdot(sau, sal)
-2

It should be:

julia> dot(sau, sal)
0

That's exactly the same issue that I could not resolve earlier...
Do you have any solution for this?

PS: I have also tested your second suggestion dotprod += dot(symmetric(A.data[j, j], sym_uplo(A.uplo), symmetric(B.data[j, j], sym_uplo(B.uplo)), which leads to the same incorrect results.

[dkarrasch]

There are transposes missing in the mixed cases, please double-check with the getindex code.

[goggle]

Ok, I think I could solve all the remaining issues, thanks @dkarrasch.
Is there anything else?

[dkarrasch]

LGTM. How about adding the Hermitian-Hermitian case also? You would need to replace all transposes by adjoints, and maybe add the real in the 2dot(...) lines. You could put all your tests in a loop over transform in (Symmetric, Hermitian), and replace all instances of Symmetric by transform, so adopting the tests should be easy.

[goggle]

I have added a second method for the Hermitian-Hermitian case and have added tests for the Hermitian case as well.

[dkarrasch]

Just a trivial comment, otherwise looks good to me. Could you maybe post an update (or a reminder) of the performance improvements by this PR?

[chriscoey]

Thanks for being so thorough on this

[goggle]

Here are some final benchmarks.

using BenchmarkTools
using LinearAlgebra
n = 5000;

Au = Symmetric(randn(n,n));
Al = Symmetric(randn(n,n), :L);
Bu = Symmetric(randn(n,n));
Bl = Symmetric(randn(n,n), :L);
Hu = Hermitian(randn(n,n) + randn(n,n)*im);
Hl = Hermitian(randn(n,n) + randn(n,n)*im, :L);
Il = Hermitian(randn(n,n) + randn(n,n)*im, :L);
Iu = Hermitian(randn(n,n) + randn(n,n)*im);

Command	Julia 1.1	Commit 8d75d17
dot(Au, Bu)	318 ms	18 ms
dot(Au, Bl)	278 ms	114 ms
dot(Al, Bu)	308 ms	112 ms
dot(Al, Bl)	388 ms	18 ms
dot(Hu, Iu)	388 ms	34 ms
dot(Hu, Il)	418 ms	135 ms
dot(Hl, Iu)	383 ms	135 ms
dot(Hl, Il)	385 ms	35 ms

@kshyatt Maybe you could add the performance label to this PR. Thanks!

[goggle]

Any news on this one? Does it need a second review? How can I proceed?

[dkarrasch]

Like a good piece of meat, this needs to hang a little... 😉 until someone with merge rights finds some time to review, and merges.

[dkarrasch]

@goggle Could you please rebase onto master/resolve the merge conflict, and then I'll merge!

[goggle]

I have rebased the branch to resolve the conflicts.
@dkarrasch Let me know if I should also squash the commits.

[dkarrasch]

Don't worry about squashing, unless this is more convenient for you while developing. Squashing will be done upon merging at the end anyway. I'll wait for tests to pass...

[dkarrasch]

Seems like one end got lost in the merge process. Add one after line 378, and then I expect it to work just fine.

[lkapelevich]

Thanks @goggle!

[goggle]

Also thanks to @dkarrasch for the great support!