# Fastest way to compute adjoint(x)\*A\*x in CUDA?

**URL:** <https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160>\
**Category:** GPU\
**Tags:** question, cuda\
**Created:** [November 2, 2024, 3:48am UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160 "2024-11-02T03:48:51Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![joseph-t-lee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joseph-t-lee/32/212353_2.png) [@joseph-t-lee](https://discourse.julialang.org/u/joseph-t-lee)\
**Post date:** [November 2, 2024, 3:48am UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/1 "2024-11-02T03:48:51Z")

</div>

I am trying to compute an expectation value of an observable, which just comes down to computing: adjoint(x) \* A \* x, where x is an N dimendional vector with ComplexF32 entries and A in an N by N matrix.

Initially A and x are both matrices. I have been converting these to CUDA arrays via

```julia
N=1024
x_d = CuArray{Float32}(undef, N)
A_d = CuArray{Float32}(undef, N, N)
copyto!(x_d, x) 
copyto!(A_d, A)

```

What would be the fastest way to compute adjoint(x) \* A \* x?

I have tried

```julia
adjoint(x_d) * A_d * x_d

```

This seems to be the slowest of all my options, as I get the impression that multiple matrix multiplications are slower. I also found that adjoint(x\_d) has a large overhead cost.

I have settled for:

```julia
tmp = A_d * x_d
Exp = dot(x_d, tmp)

```

Is there a faster option to do this? Also, is there any way to use sparse arrays? And is it necessarily faster?

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/224680_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [November 2, 2024, 6:32am UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/2 "2024-11-02T06:32:38Z")

</div>

There’s a 3-argument `dot(y, A, x)` function in `LinearAlgebra`, you could try it and see if it is overloaded for `CuArray`s?

---

<div class="post-metadata">

**Author:** ![joseph-t-lee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joseph-t-lee/32/212353_2.png) [@joseph-t-lee](https://discourse.julialang.org/u/joseph-t-lee)\
**Post date:** [November 2, 2024, 5:53pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/3 "2024-11-02T17:53:46Z")

</div>

It looks like dot(y, A, x) has a huge overhead cost and is quite slow. I was comparing the following:

```julia
d = 1048
x = rand(ComplexF32, d)
A = rand(ComplexF32, d, d)
x_d = CuArray{ComplexF32}(undef, d)
A_d = CuArray{ComplexF32}(undef, d, d)
copyto!(x_d, x)
copyto!(A_d, A)

@time tmp = A_d * x_d
@time res1 = dot(x_d, tmp)
@time res2 = adjoint(x_d) * A_d * x_d
@time res3 = dot(x_d, A_d, x_d)
println()
@btime tmp = $A_d * $x_d
@btime res1 = dot($x_d, $tmp)
@btime res2 = adjoint($x_d) * $A_d * $x_d
@btime res3 = dot($x_d, $A_d, $x_d)

```

yielding:

```julia
0.000144 seconds (23 allocations: 496 bytes)
 0.001023 seconds (15 allocations: 240 bytes)
 0.000358 seconds (39 allocations: 752 bytes)
26.118309 seconds (15.39 M allocations: 301.947 MiB, 0.15% gc time)

 5.660 μs (23 allocations: 496 bytes)
 20.800 μs (14 allocations: 224 bytes)
 87.400 μs (37 allocations: 720 bytes)
 26.185 s (15390956 allocations: 301.95 MiB)

```

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [November 2, 2024, 5:59pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/4 "2024-11-02T17:59:44Z")

</div>

I don’t think three-arg `dot` is implemented for `CuArray`s (at least it didn’t work for me), so you’re probably falling back to a generic loop with scalar indexing, hence the massive overhead.

Try adding some parentheses to reassociate, i.e., try `x' * (A * x)`. That’s the fastest version on my end, unless you go down the route of preallocating a temp vector.

Btw. note that

```julia
tmp = A * x
dot(x, tmp)

```

can be simplified to the equivalent `dot(x, A * x)`.

---

<div class="post-metadata">

**Author:** ![joseph-t-lee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joseph-t-lee/32/212353_2.png) [@joseph-t-lee](https://discourse.julialang.org/u/joseph-t-lee)\
**Post date:** [November 2, 2024, 6:05pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/5 "2024-11-02T18:05:38Z")

</div>

I’m still finding that

```julia
tmp = A * x
dot(x, tmp)

```

computes in the same order of magnitude as

```julia
dot(x, A * x)

```

with @time, and the former is faster with @btime.

I don’t know if I fully understand how GPU computing works. When I write a line like ` y * (A * x)` does my computer interpret this as do A\*x first, then y \* (Ax)?

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [November 2, 2024, 6:57pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/6 "2024-11-02T18:57:17Z")

</div>

> [@joseph-t-lee](#):
>
> I’m still finding that
> 
> ```julia
> tmp = A * x
> dot(x, tmp)
> 
> ```
> 
> computes in the same order of magnitude as
> 
> ```julia
> dot(x, A * x)
> 
> ```
> 
> with @time, and the former is faster with @btime.

That’s strange. Here’s what I get:

```julia-repl
julia> using CUDA, LinearAlgebra, BenchmarkTools

julia> N = 1024;

julia> x = CUDA.rand(ComplexF32, N);

julia> A = CUDA.rand(ComplexF32, N, N);

julia> @btime begin
           tmp = $A * $x
           dot($x, tmp)
       end;
  118.237 μs (37 allocations: 720 bytes)

julia> @btime dot($x, $A * $x);
  116.874 μs (37 allocations: 720 bytes)

```

> [@joseph-t-lee](#):
>
> When I write a line like ` y * (A * x)` does my computer interpret this as do A\*x first, then y \* (Ax)?

Yes. This is just what parentheses do, regardless of whether `A` or `x` are GPU arrays.

---

<div class="post-metadata">

**Author:** ![RomeoV](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romeov/32/37687_2.png) [@RomeoV](https://discourse.julialang.org/u/RomeoV)\
**Post date:** [November 2, 2024, 7:20pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/7 "2024-11-02T19:20:23Z")

</div>

FYI I don’t think CUDA itself supports x^T A x directly. So regardless how you implement it in Julia, in the end you’ll have to compute something like x^T (A x) anyways.

Also x^T y will always be “slow” on GPUs as you’re reducing a vector, which GPUs can’t do particularly well. If you need more speedup, you may look into whether you can batch multiple of those products together (if your program works that way).  
I’m not sure at what size it’s worth it to transfer to CPU and reduce there…

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [November 2, 2024, 7:27pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/8 "2024-11-02T19:27:07Z")

</div>

> [@RomeoV](#):
>
> x^T y will always be “slow” on GPUs as you’re reducing a vector, which GPUs can’t do particularly well.

Is this true? A GPU can do parallel tree reduction with complexity O(log N), while on a CPU this is always an O(N) operation (tree reduction can still reduce accumulated roundoff error, but not complexity). Sure you probably need a pretty large N to overcome the slower clock speed of the GPU, but it seems to have the advantage of scaling. (Caveat: haven’t run any benchmarks.)

---

<div class="post-metadata">

**Author:** ![joseph-t-lee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joseph-t-lee/32/212353_2.png) [@joseph-t-lee](https://discourse.julialang.org/u/joseph-t-lee)\
**Post date:** [November 2, 2024, 7:34pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/9 "2024-11-02T19:34:21Z")

</div>

I redid the benchmarking - doing the multiplication in 2 steps is equally as fast as dot($x, $A \* $x), so I messed something up beforehand.

I was also under the impression that any type of matrix multiplication, dot product included, was faster on GPU. Though admittedly I am quite new to GPU computing and have no idea what the details are.

---

<div class="post-metadata">

**Author:** ![RomeoV](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romeov/32/37687_2.png) [@RomeoV](https://discourse.julialang.org/u/RomeoV)\
**Post date:** [November 2, 2024, 8:04pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/10 "2024-11-02T20:04:09Z")

</div>

> Is this true?

Let me be more clear. Reducing on the GPU is slow compared to other operations on the GPU. Let’s compare summing a large matrix along one dimension to summing (reducing) a vector on a GPU.  
~~Even though summing the matrix takes a lot more mathematical operations, it is faster than reducing the vector. The reason for this is what you mentioned.~~ redacted as per below.

```julia
using CUDA
x = cu(rand(10_000))
A = cu(rand(10_000, 10_000))

# EDIT: These numbers seem to be misleading.
julia> CUDA.@profile sum(x)
Host-side activity: calling CUDA APIs took 80.11 µs (38.44% of the trace)
# 80.11 µs
[...]
julia> CUDA.@profile sum(A; dims=1)
Host-side activity: calling CUDA APIs took 36.48 µs (0.79% of the trace)
# 36.48 µs
[...]

```

So the reduction is “slow” in the sense that it doesn’t compare favorably to other operations on the GPU.

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [November 2, 2024, 8:21pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/11 "2024-11-02T20:21:47Z")

</div>

Not sure this is how to read the output of `CUDA.@profile`, unless I’m misunderstanding your claim. Here’s what plain timings show:

```julia-repl
julia> x = CUDA.rand(10000);

julia> A = CUDA.rand(10000, 10000);

julia> @btime CUDA.@sync sum($x);
  91.865 μs (95 allocations: 2.88 KiB)

julia> @btime CUDA.@sync sum($A; dims=1);
  110.844 ms (34 allocations: 976 bytes)

julia> @btime CUDA.@sync sum($A; dims=2);
  5.980 ms (34 allocations: 976 bytes)

```

---

<div class="post-metadata">

**Author:** ![RomeoV](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romeov/32/37687_2.png) [@RomeoV](https://discourse.julialang.org/u/RomeoV)\
**Post date:** [November 2, 2024, 8:35pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/12 "2024-11-02T20:35:02Z")

</div>

Yes you’re right, I am getting the same as you when measuring in another way. Thanks for the correction. I suppose CUDA.@profile doesn’t report what I thought, and also I had a thinking error.  
When reducing a NxN matrix the GPU needs to reduce N times. When reducing a vector only log2(N) times. So we would expect `sum(x)` to be faster.

However, it still stays true that 1D reductions on the GPU are “slow” in the sense that GPUs can’t exploit their full parallelism.

I’ve also tried to benchmark `x' * A * x` vs `collect(x)' * collect(A * x)` where `A` and `x` are GPU arrays. However, I’m not sure if the results I’m getting are true. (My benchmark indicates that they take about the same time, across different sizes).

```julia
N = 10_000
x=cu(rand(Float32, N)); A=cu(rand(Float32, N, N))
julia> @btime CUDA.@sync x' * (A*x);
  2.379 ms (39 allocations: 752 bytes)

julia> @btime CUDA.@sync collect(x)' * collect(A*x);
  2.435 ms (44 allocations: 157.12 KiB)

```

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [November 2, 2024, 8:39pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/13 "2024-11-02T20:39:15Z")

</div>

> [@RomeoV](#):
>
> However, it still stays true that 1D reductions on the GPU are “slow” in the sense that GPUs can’t exploit their full parallelism.

At `N = 10000`, my dinky laptop GPU certainly isn’t of any help compared to the CPU. Continuing from the above:

```julia-repl
julia> x_h = Array(x);

julia> @btime sum($x_h); # Compare GPU: 91.865 μs
  1.364 μs (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [November 2, 2024, 8:41pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/14 "2024-11-02T20:41:18Z")

</div>

> [@RomeoV](#):
>
> My benchmark indicates that they take about the same time, across different sizes

Probably because the timings are completely dominated by the matrix-vector product at every size you’re checking?

---

<div class="post-metadata">

**Author:** ![RomeoV](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/romeov/32/37687_2.png) [@RomeoV](https://discourse.julialang.org/u/RomeoV)\
**Post date:** [November 2, 2024, 8:44pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/15 "2024-11-02T20:44:35Z")

</div>

Yeah that seems about right.

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [November 2, 2024, 8:45pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/16 "2024-11-02T20:45:11Z")

</div>

> [@RomeoV](#):
>
> When reducing a NxN matrix the GPU needs to reduce N times.

For an N x N matrix you can do N tree reductions in parallel (or some interpolation between the two extremes). I don’t know what `sum` is actually doing here, but it may be limited more by the number of threads that can actually run in parallel, rather than algorithmic limitations, depending on hardware.

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [November 2, 2024, 8:54pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/17 "2024-11-02T20:54:52Z")

</div>

> [@joseph-t-lee](#):
>
> Also, is there any way to use sparse arrays? And is it necessarily faster?

I wanted to circle back to this.

1. Sparse arrays are generally only useful if they are _really_ sparse, which typically means that only a few percent of the entries are nonzero. If your problem can be instantiated at multiple sizes N, the number of nonzero elements in the matrix should probably scale slower than N^2.
2. You can use sparse arrays both on CPU and GPU. I’ve never tried it on GPU, but the relevant Cuda module is `CUDA.CUSPARSE`. However, I wouldn’t expect to see a speedup of GPU compared to CPU as often for sparse problems. GPUs really excel at large, dense matrix-matrix multiplication—the more elements and the more work per element, the better.

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [November 2, 2024, 9:06pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/18 "2024-11-02T21:06:51Z")

</div>

Getting back to the actual problem: If you’re calculating lots of expectation values at the same system size, you can save a some time by preallocating a temporary array to hold `A * x`, so you don’t have to allocate a new one every time. It’s not a huge speedup, and relatively speaking it gets smaller as the system size increases (since `A * x` scales as N^2 while the allocation scales as N), but it makes a measurable difference at N = 1024 on my machine:

```julia-repl
julia> using CUDA, LinearAlgebra, BenchmarkTools

julia> x = CUDA.rand(ComplexF32, 1024);

julia> A = CUDA.rand(ComplexF32, 1024, 1024);

julia> tmp = similar(x);

julia> @btime dot($x, $A * $x); # not using tmp
  118.572 μs (37 allocations: 720 bytes)

julia> @btime dot($x, mul!($tmp, $A, $x)); # using tmp
  111.612 μs (32 allocations: 512 bytes)

```

---

<div class="post-metadata">

**Author:** ![joseph-t-lee](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joseph-t-lee/32/212353_2.png) [@joseph-t-lee](https://discourse.julialang.org/u/joseph-t-lee)\
**Post date:** [November 2, 2024, 9:10pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/19 "2024-11-02T21:10:48Z")

</div>

Thanks for the input everyone. So my takeaways are the following:

- Matrix multiplication is faster on the GPU (of course for small N it is probably not worthwhile). I can be smart about my parenthesis etc. to tell my computer what order to do things in. (I am assuming it is true that A\*x is also still faster even if x is a vector? Though maybe I am misinterpreting the benchmarking above).
- I’m not sure if I fully understood the discussion above, but I take it that sums over vectors (i.e. x^T A x type calculations) won’t be that beneficial on a GPU.

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [November 2, 2024, 9:28pm UTC](https://discourse.julialang.org/t/fastest-way-to-compute-adjoint-x-a-x-in-cuda/122160/20 "2024-11-02T21:28:54Z")

</div>

> [@joseph-t-lee](#):
>
> I am assuming it is true that A\*x is also still faster even if x is a vector?

It can definitely be. You should verify with benchmarks on the machine you’re using, with something that looks like your problem, but here’s my laptop with a dense random matrix (remember to use `CUDA.@sync`!):

```julia-repl
julia> using CUDA, LinearAlgebra, BenchmarkTools

julia> x = CUDA.rand(ComplexF32, 1024);

julia> A = CUDA.rand(ComplexF32, 1024, 1024);

julia> tmp = similar(x);

julia> x_h, A_h, tmp_h = Array(x), Array(A), Array(tmp); # copy to CPU

julia> @btime CUDA.@sync mul!($tmp, $A, $x); # GPU
  101.661 μs (18 allocations: 288 bytes)

julia> @btime mul!($tmp_h, $A_h, $x_h); # CPU
  236.681 μs (0 allocations: 0 bytes)

```

> [@joseph-t-lee](#):
>
> I take it that sums over vectors (i.e. x^T A x type calculations) won’t be that beneficial on a GPU.

This is specifically about linear reductions such as the dot product between two vectors, `x' * tmp` or equivalently `dot(x, tmp)`, which might be slower on the GPU than the CPU.

```julia-repl
julia> @btime CUDA.@sync dot($x, $tmp); # GPU
  23.040 μs (14 allocations: 224 bytes)

julia> @btime dot($x_h, $tmp_h); # CPU
  307.195 ns (0 allocations: 0 bytes)

```

But for large N, the expectation value x^T A x is completely dominated by the matrix-vector product `A * x`. So if the matrix-vector product gets a significant speedup on GPU, it’s worth using the GPU even if the final step is a little slower there.

```julia-repl
julia> @btime CUDA.@sync dot($x, mul!($tmp, $A, $x)); # GPU
  112.486 μs (32 allocations: 512 bytes)

julia> @btime dot($x_h, mul!($tmp_h, $A_h, $x_h)); # CPU
  244.372 μs (0 allocations: 0 bytes)

```

You can also try a best-of-both-worlds strategy where you do A \* x on the GPU and then copy the result to the CPU to compute the dot product, but on my machine this is slower due to the overhead of moving data from the GPU to the CPU:

```julia-repl
julia> @btime begin
           mul!($tmp, $A, $x) # GPU
           dot(Array($x), Array($tmp)) # CPU
       end;
  131.419 μs (30 allocations: 16.69 KiB)

```