# Efficient and in-place computation of A + adjoint(A)
**URL:** https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290
**Category:** Performance
**Tags:** question
**Created:** [January 21, 2023, 12:24am UTC](https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290 "2023-01-21T00:24:13Z")
**Posts on this page:** 15
**Page:** 1
### Author:  [@Ronan\_Gautier](https://discourse.julialang.org/u/Ronan_Gautier)
#### Post date: [January 21, 2023, 12:24am UTC](https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290/1 "2023-01-21T00:24:13Z")
If `A` is any complex-valued matrix, is there a more efficient way of computing `A + adjoint(A)` in-place than the following ?
```julia
function hermitian(A, A_tmp)
adjoint!(A_tmp, A)
A .+= A_tmp
end
```
Benchmarking with
```julia
using LinearAlgebra
using BenchmarkTools
const N = 200
A = rand(ComplexF64, N, N)
A_tmp = zeros(ComplexF64, N, N)
@btime hermitian(A, A_tmp)
```
yields `70.243 μs (0 allocations: 0 bytes)`, which seems quite slow for such a small matrix. Is there any other more efficient implementation ?
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [January 21, 2023, 12:36am UTC](https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290/2 "2023-01-21T00:36:33Z")
See [https://github.com/JuliaLang/julia/pull/31836](https://github.com/JuliaLang/julia/pull/31836)
---
### Author:  [@Ronan\_Gautier](https://discourse.julialang.org/u/Ronan_Gautier)
#### Post date: [January 21, 2023, 1:17am UTC](https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290/3 "2023-01-21T01:17:12Z")
Thanks for the link. So, just for anyone else who might be interested. It seems that the following is the most efficient if you want the full matrix to be updated at the end.
```julia
function hermitian!(A)
@inbounds for j in 1:N
A[j, j] = 2 * real(A[j, j])
for i in 1:(j-1)
A[i, j] = A[i, j] + conj(A[j, i])
end
end
@inbounds for j in 1:N
for i in j+1:N
A[i, j] = conj(A[j, i])
end
end
end
```
If you’re only interested in the upper diagonal part of the matrix, you can skip the second for loop, and gain a factor of 2 on efficiency.
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [January 21, 2023, 1:43am UTC](https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290/4 "2023-01-21T01:43:14Z")
A little bit faster (on my machine) is:
```julia
function hermitian!(A::Matrix{Complex{T}}) where T
( N = size(A,1) ) == size(A, 2) || error("Matrix must be square")
@inbounds for j in 1:N
A[j, j] = 2 * real(A[j, j])
for i in 1:(j-1)
r1, i1 = reim(A[i,j])
r2, i2 = reim(A[j,i])
rc, ic = r1+r2, i1-i2
A[i, j], A[j, i] = Complex{T}(rc,ic), Complex{T}(rc,-ic)
end
end
end
```
Once `A[j,i]` is touched in the first loop, it is better to write to it while it is hot in the cache.
---
### Author:  [@Seif\_Shebl](https://discourse.julialang.org/u/Seif_Shebl)
#### Post date: [January 22, 2023, 6:19pm UTC](https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290/5 "2023-01-22T18:19:39Z")
This simpler version would be as fast as Dan’s solution:
```julia
function hermitian!(A)
(N = size(A,1)) == size(A,2) || error("Matrix must be square")
@inbounds for i = 1:N, j = i:N
A[i, j] = A[i, j] + A[j, i]'
A[j, i] = A[i, j]'
end
end
```
---
### Author:  [@uniment](https://discourse.julialang.org/u/uniment)
#### Post date: [January 23, 2023, 12:43am UTC](https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290/6 "2023-01-23T00:43:11Z")
Note that if `A` is an `OffsetArray`, Julia may terminate with this error:
```plaintext
Please submit a bug report with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks.
Exception: EXCEPTION_ACCESS_VIOLATION at 0x7ffba590f6c1 -- RtlAllocateHeap at C:\WINDOWS\SYSTEM32\ntdll.dll (unknown line)
in expression starting at none:0
RtlAllocateHeap at C:\WINDOWS\SYSTEM32\ntdll.dll (unknown line)
RtlAllocateHeap at C:\WINDOWS\SYSTEM32\ntdll.dll (unknown line)
malloc at C:\WINDOWS\System32\msvcrt.dll (unknown line)
aligned_offset_malloc at C:\WINDOWS\System32\msvcrt.dll (unknown line)
```
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [January 23, 2023, 7:57am UTC](https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290/7 "2023-01-23T07:57:47Z")
> [@Seif\_Shebl](#):
>
> `@inbounds for i = 1:N, j = i:N`
`@inbounds` together with 1-based indexing is not advisable. Use `axes`. Many have argued lately that it’s fine to use `1:size()`, but then remove `@inbounds`.
---
### Author:  [@photor](https://discourse.julialang.org/u/photor)
#### Post date: [January 23, 2023, 8:55am UTC](https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290/8 "2023-01-23T08:55:22Z")
```julia
for i in axes(A,1)
```
but how to index j?
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [January 23, 2023, 9:34am UTC](https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290/9 "2023-01-23T09:34:01Z")
```julia
for i in axes(A, 1), j in i:lastindex(A, 2)
```
Or, simply, the way it is already, but _without_ the `@inbounds` annotation.
`@inbounds` is for when the indices are _provably_ inbounds. That could be achieved by using `eachindex`, `axes`, etc. Or by only accepting `Array` as input type, or by calling `Base.require_one_based_indexing` if you strongly prefer `1:N` indexing.
---
### Author:  [@uniment](https://discourse.julialang.org/u/uniment)
#### Post date: [January 24, 2023, 2:45am UTC](https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290/10 "2023-01-24T02:45:40Z")
Somebody came up with a really clever solution to this problem, it just needs to be fleshed out and implemented
> [@Is OffsetArrays.jl a poison pill?](https://discourse.julialang.org/t/is-offsetarrays-jl-a-poison-pill/85188/21):
>
> But maybe there could be some very simple and universal way to get a one-based view into any `AbstractArray`? `OffsetArray` has `no_offset_view`, but should all `AbstractArray`s have a method like that?
---
### Author:  [@jishnub](https://discourse.julialang.org/u/jishnub)
#### Post date: [January 24, 2023, 5:56am UTC](https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290/11 "2023-01-24T05:56:00Z")
Without restarting the whole offset array discussion again, `reshape(A, size(A))` should return a 1-based view for all arrays
```julia
julia> A = zeros(3:4, 4:5)
2×2 OffsetArray(::Matrix{Float64}, 3:4, 4:5) with eltype Float64 with indices 3:4×4:5:
0.0 0.0
0.0 0.0
julia> reshape(A, size(A))
2×2 Matrix{Float64}:
0.0 0.0
0.0 0.0
```
There might be missing cases that need to be fixed, but this is what the function is expected to return.
---
### Author:  [@uniment](https://discourse.julialang.org/u/uniment)
#### Post date: [January 24, 2023, 7:55am UTC](https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290/12 "2023-01-24T07:55:00Z")
For those of us who like 1-based indexing (🙋♂️) but want performance and genericism without ripping up our code and making it unrecognizeable, shall this be the preferred idiom?
```julia
hermitian!(A::AbstractArray) = let A=reshape(A, size(A))
(N = size(A,1)) == size(A,2) || error("Matrix must be square")
@inbounds for i = 1:N, j = i:N
A[i, j] = A[i, j] + A[j, i]'
A[j, i] = A[i, j]'
end
end
```
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [January 24, 2023, 9:51am UTC](https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290/13 "2023-01-24T09:51:11Z")
> [@jishnub](#):
>
> `A = zeros(3:4, 4:5)`
Is this the notation for creating an OffsetArray of zeros? Isn’t that type piracy?
---
### Author:  [@uniment](https://discourse.julialang.org/u/uniment)
#### Post date: [January 24, 2023, 10:04am UTC](https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290/14 "2023-01-24T10:04:54Z")
I argue that `OffsetArray`s should be added to `Base`, in which case it wouldn’t be piracy anymore 😉
---
### Author:  [@jishnub](https://discourse.julialang.org/u/jishnub)
#### Post date: [January 24, 2023, 10:14am UTC](https://discourse.julialang.org/t/efficient-and-in-place-computation-of-a-adjoint-a/93290/15 "2023-01-24T10:14:04Z")
Yes, it’s type-piracy, and is [expected to be removed](https://github.com/JuliaArrays/OffsetArrays.jl/issues/306) in the next breaking release of OffsetArrays.