# Element wise multiplication of lower triangular matrices
**URL:** https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252
**Category:** New to Julia
**Tags:** linearalgebra
**Created:** [August 27, 2023, 8:36am UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252 "2023-08-27T08:36:22Z")
**Posts on this page:** 18
**Page:** 1
### Author:  [@ahmed](https://discourse.julialang.org/u/ahmed)
#### Post date: [August 27, 2023, 8:36am UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/1 "2023-08-27T08:36:22Z")
I’m confused by how the element wise multiplication of lower triangular matrices is not faster than the element wise multiplication of the completely full matrices. This behavior is shown in the following MWE
```julia
using LinearAlgebra, BenchmarkTools
A = rand(ComplexF64, 1000, 1000);
B = rand(ComplexF64, 1000, 1000);
function f(A, B)
UnitLowerTriangular(A) .* UnitLowerTriangular(B)
end
function f2(A, B)
return A .* B
end
```
The benchmark of the two functions give
```julia
@btime f(A, B)
4.017 ms (2 allocations: 15.26 MiB)
```
```julia
@btime f2(A, B)
3.719 ms (2 allocations: 15.26 MiB)
```
The number of operations in the lower triangular matrices should be close to half of that in the complete matrices; so I’d think it should take half the time?
---
### Author:  [@Mason](https://discourse.julialang.org/u/Mason)
#### Post date: [August 27, 2023, 8:48am UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/2 "2023-08-27T08:48:02Z")
The problem is that for _most_ operations involving symmetric or triangular matrices, yes you can get a 2x speedup by only doing half the operations, but it would come at the cost of losing the typically ~8x speedup you can get from SIMD.
For this reason, triangular and symmetric matrices are typically actually stored just like regular matrices, but the lower or upper triangle are treated as ‘junk’ bits which are ignored.
So if we used denser linear storage yes we could get 2x faster element-wise multiplication, but things like matrix-multiplication would suffer immensely.
---
### Author:  [@ahmed](https://discourse.julialang.org/u/ahmed)
#### Post date: [August 27, 2023, 9:01am UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/3 "2023-08-27T09:01:10Z")
Thanks @Mason for your insight. Actually the reason I’m asking this question is because I’m doing an element wise multiplication of two Hermitian matrices and since the answer will be another Hermitian matrix so I was hoping I could save almost half the time this way.
I’m still not really clear on why this speed up can’t be done.
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [August 27, 2023, 11:13am UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/4 "2023-08-27T11:13:03Z")
Maybe one could build a data structure with both representations.
One will be used for matrix operations while the other for element wise operations.
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [August 27, 2023, 1:43pm UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/5 "2023-08-27T13:43:09Z")
> [@Mason](#):
>
> it would come at the cost of losing the typically ~8x speedup you can get from SIMD.
There’s no reason you can’t have it both ways, at least if you help the simd manually. I presume there’s a challenge getting it to happen automatically, though.
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [August 27, 2023, 2:54pm UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/6 "2023-08-27T14:54:24Z")
I think @Mason meant the loading of the data is not from a contiguous memory.
Hence the performance benefit is much lower as you need to do multiple non `SIMD` loads into a vector and then do the `SIMD` processing (Multiplication).
---
### Author:  [@dmbates](https://discourse.julialang.org/u/dmbates)
#### Post date: [August 27, 2023, 3:00pm UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/7 "2023-08-27T15:00:46Z")
I think the point is that matrix multiplication is one of the most highly optimized operations in any custom BLAS suite. At a minimum Julia uses OpenBlas but there are other options such as [MKL](https://github.com/JuliaLinearAlgebra/MKL.jl) for Intel processors and [AppleAccelerate](https://github.com/JuliaLinearAlgebra/AppleAccelerate.jl) for Apple M-series processors. A lot of effort went into those implementations by people at the respective companies working closely with the hardware folks. Reproducing such intensely optimized code for triangular or Hermitian matrices would be a huge undertaking.
There is the [RectangularFullPacked](https://github.com/JuliaLinearAlgebra/RectangularFullPacked.jl) package available. It uses a reduced-size representation of triangular or Hermitian matrices that doesn’t lose too much performance relative to having a full size square array as a backing store.
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [August 27, 2023, 5:58pm UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/8 "2023-08-27T17:58:10Z")
> [@RoyiAvital](#):
>
> Hence the performance benefit is much lower as you need to do multiple non `SIMD` loads into a vector and then do the `SIMD` processing (Multiplication).
You can do all simd loads and process everything with simd. You must just accept that not each load will be contiguous with the previous one, and that there will be a slight redundancy in how much data you must process.
> [@dmbates](#):
>
> I think the point is that matrix multiplication is one of the most highly optimized operations in any custom BLAS suite
This particular case was elementwise multiplication, though, not matrix product.
---
### Author:  [@Mason](https://discourse.julialang.org/u/Mason)
#### Post date: [August 27, 2023, 7:04pm UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/9 "2023-08-27T19:04:55Z")
> [@DNF](#):
>
> There’s no reason you can’t have it both ways, at least if you help the simd manually. I presume there’s a challenge getting it to happen automatically, though.
I never meant that if was impossible to get this 2x speedup for triangular structures (in fact, I did say that it was possible for element wise stuff). My point was just that the algorithms and data structures are all designed in a format that makes getting this last 2x kinda inconvenient and awkward.
Someone could definitely do it, it’s just that elementwise operations arent really something people typically care enough about to optimize that hard.
---
### Author:  [@dmbates](https://discourse.julialang.org/u/dmbates)
#### Post date: [August 27, 2023, 7:16pm UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/10 "2023-08-27T19:16:02Z")
You’re right. This is what happens when I reply to messages before coffee.
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [August 29, 2023, 7:21am UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/11 "2023-08-29T07:21:49Z")
> [@DNF](#):
>
> You can do all simd loads and process everything with simd. You must just accept that not each load will be contiguous with the previous one, and that there will be a slight redundancy in how much data you must process.
It is not only contiguous, it will be not aligned.
Unless you do indices tricks in the loop which then means you will need to do everything in hand as the compiler will see it as data dependency.
So probably you have to create a different data structure optimized for this.
Which will cause issues in the cases of element wise multiplication with other structures and other operations.
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [August 29, 2023, 7:44am UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/12 "2023-08-29T07:44:53Z")
> [@RoyiAvital](#):
>
> It is not only contiguous, it will be not aligned.
Not sure what you mean. If the array is aligned and you have regular numbers (e.g. 32 or 64-bit) you can start your load at any arbitrary position in the array, except you must be careful near the end (I’ve been spending the past few weeks jumping around arrays with avx loads exactly like that).
> [@RoyiAvital](#):
>
> Unless you do indices tricks in the loop which then means you will need to do everything in hand as the compiler will see it as data dependency.
You need to calculate when and how far to skip for each column. Depending on how hard this is, it may end up costing more than what you save, but I think it should work well for large matrices.
As @Mason said, the benefits may not be great enough to justify the extra work and complexity, in the general case.
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [August 29, 2023, 12:37pm UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/13 "2023-08-29T12:37:42Z")
Just for fun, I implemented my idea using SIMD.jl. It only works for real-valued matrices, complex values are a different kettle of fish.
For 1000x1000 matrices, it is approximately 2x as fast as a regular elementwise multiplication over the (dense) parent array, pretty close to what I predicted. (I also tried LoopVectorization.jl, which offered no speedup.) It is not fully optimal, nor is it really ‘safe’ with respect to all sizes and indexing schemes (it will probably fail for very small matrices, for example):
```julia
function simdmul!(C::LowerTriangular{T}, A::LowerTriangular{T}, B::LowerTriangular{T}) where {T<:Real}
size(A) == size(B) == size(C) || error("Incompatible sizes.")
N = div(512, 8 * sizeof(T))
lowermul!(parent(C), parent(A), parent(B), Val(N))
end
function lowermul!(C, A, B, ::Val{N}) where {N}
lane = VecRange{N}(firstindex(A, 1))
(J, K) = size(C)
J_ = J
skip = 0
for _ in 1:K-1
for j in skip:N:J_
lanej = lane + j
C[lanej] = A[lanej] * B[lanej]
end
skip += J + 1
J_ += J
end
C[end] = A[end] * B[end]
return C
end
```
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [August 29, 2023, 12:57pm UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/14 "2023-08-29T12:57:05Z")
Just for fun, I implemented my idea using SIMD.jl. It only works for real-valued matrices, complex values are a different kettle of fish.
For 1000x1000 matrices, it is approximately 2x as fast as a regular elementwise multiplication over the (dense) parent array, pretty close to what I predicted. I also tried LoopVectorization It is not fully optimal, nor is it really ‘safe’ with respect to all sizes and indexing schemes (it will probably fail for very small matrices, for example):
```julia
function simdmul!(C::LowerTriangular{T}, A::LowerTriangular{T}, B::LowerTriangular{T}) where {T<:Real}
size(A) == size(B) == size(C) || error("Incompatible sizes.")
N = div(512, 8 * sizeof(T))
lowermul!(parent(C), parent(A), parent(B), Val(N))
end
function lowermul!(C, A, B, ::Val{N}) where {N}
lane = VecRange{N}(firstindex(A, 1))
(J, K) = size(C)
J_ = J
skip = 0
for _ in 1:K-1
for j in skip:N:J_
lanej = lane + j
C[lanej] = A[lanej] * B[lanej]
end
skip += J + 1 # <- this is the trick that saves operations
J_ += J
end
C[end] = A[end] * B[end]
return C
end
```
**Edit:** Actually, you don’t need manual SIMD, you can just ~~ask LoopVectorization to vectorize the inner loop for you~~ do the completely simple naive thing, leading to cleaner code:
```julia
function automul!(C::LowerTriangular{T}, A::LowerTriangular{T}, B::LowerTriangular{T}) where {T<:Real}
size(A) == size(B) == size(C) || error("Incompatible sizes.")
lowermul!(parent(C), parent(A), parent(B))
return C
end
function lowermul!(C, A, B)
for k in axes(C, 2)
for j in k:lastindex(C, 1)
C[j, k] = A[j, k] * B[j, k]
end
end
return C
end
```
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [August 29, 2023, 2:24pm UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/15 "2023-08-29T14:24:50Z")
> [@DNF](#):
>
> **Edit:** Actually, you don’t need manual SIMD, you can just ~~ask LoopVectorization to vectorize the inner loop for you~~ do the completely simple naive thing, leading to cleaner code:
>
> ```julia
> function automul!(C::LowerTriangular{T}, A::LowerTriangular{T}, B::LowerTriangular{T}) where {T<:Real}
> size(A) == size(B) == size(C) || error("Incompatible sizes.")
> lowermul!(parent(C), parent(A), parent(B))
> return C
> end
>
> function lowermul!(C, A, B)
> for k in axes(C, 2)
> for j in k:lastindex(C, 1)
> C[j, k] = A[j, k] * B[j, k]
> end
> end
> return C
> end
>
> ```
I was under the impression Julia won’t generate `SIMD` code for triangle loop.
Can you verify? It seems strange.
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [August 29, 2023, 2:40pm UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/16 "2023-08-29T14:40:23Z")
Here’s the weird thing. If I do
```julia
function tmul!(C, A, B)
@turbo for i in eachindex(A, B, C)
C[i] = A[i] * B[i]
end
return C
end
```
and call it with regular 1000x1000 arrays, `@code_llvm` shows plenty of `load <4 x double>`. The vanilla triangular loop shows _no_ vectorization in `@code_llvm`. And _still_, it is _twice_ as fast.
What’s going on?
---
### Author:  [@mikmoore](https://discourse.julialang.org/u/mikmoore)
#### Post date: [August 29, 2023, 6:06pm UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/17 "2023-08-29T18:06:42Z")
> [@DNF](#):
>
> And _still_, it is _twice_ as fast.
>
> What’s going on?
Is this a symptom of memory-bound operations? SIMD doesn’t make memory any faster while the triangular loop only pulls ~half the memory.
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [August 29, 2023, 6:20pm UTC](https://discourse.julialang.org/t/element-wise-multiplication-of-lower-triangular-matrices/103252/18 "2023-08-29T18:20:32Z")
Yeah, that’s it. Reducing the size to 20x20 changes the picture.