# Absolute value of a matrix

**URL:** <https://discourse.julialang.org/t/absolute-value-of-a-matrix/98301>\
**Category:** General Usage\
**Tags:** performance, linearalgebra\
**Created:** [May 4, 2023, 7:10am UTC](https://discourse.julialang.org/t/absolute-value-of-a-matrix/98301 "2023-05-04T07:10:23Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![e3c6](https://avatars.discourse-cdn.com/v4/letter/e/e79b87/32.png) [@e3c6](https://discourse.julialang.org/u/e3c6)\
**Post date:** [May 4, 2023, 7:10am UTC](https://discourse.julialang.org/t/absolute-value-of-a-matrix/98301/1 "2023-05-04T07:10:23Z")

</div>

Let `A` be a symmetric (or Hermitian) matrix. I want to compute `abs(A)`. By this I mean taking the absolute value of the eigenvalues. See for instance [linear algebra - Understanding the absolute value of a matrix. - Mathematics Stack Exchange](https://math.stackexchange.com/questions/2829437/understanding-the-absolute-value-of-a-matrix).

I was surprised that `Base.abs(A)` does not work. I implemented a few alternatives here:

```julia
absmat1(A::AbstractMatrix) = sqrt(A'A)
absmat2(A::AbstractMatrix) = sqrt(Hermitian(A'A))
function absmat3(A::AbstractMatrix)
    F = eigen(A)
    return F.vectors * Diagonal(abs.(F.values)) * F.vectors'
end

```

Some benchmarks (in the style of [Fast `diag(A' * B * A)` - #6 by Mason](https://discourse.julialang.org/t/fast-diag-a-b-a/98216/6)):

```julia
let n = 1000
    _A = randn(n, n)
    A = _A + _A'
    B = Hermitian(A)
    for f ∈ (absmat1, absmat2, absmat3)
        print(f, " ")
        @btime $f($A)
        print(f, " (sym) ")
        @btime $f($B)
    end
end;
# absmat1 133.304 ms (21 allocations: 46.14 MiB)
# absmat1 (sym) 139.904 ms (23 allocations: 53.77 MiB)
# absmat2 132.386 ms (21 allocations: 46.14 MiB)
# absmat2 (sym) 139.972 ms (23 allocations: 53.77 MiB)
# absmat3 147.802 ms (20 allocations: 38.51 MiB)
# absmat3 (sym) 147.689 ms (18 allocations: 38.51 MiB)

```

Is there a faster way?

I would also suggest `abs(A)` be added to `Base`.

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [May 4, 2023, 8:41am UTC](https://discourse.julialang.org/t/absolute-value-of-a-matrix/98301/2 "2023-05-04T08:41:18Z")

</div>

I would look into SVD functions in `LinearAlgebra`

> **[Singular value](https://en.wikipedia.org/wiki/Singular_value?wprov=sfti1)**
>
> In mathematics, in particular functional analysis, the singular values of a compact operator 
>   
>     
>       
> T
> :
> X
> →
> Y
>       
>     
> {\\displaystyle T:X\\rightarrow Y}
>   
> acting between Hilbert spaces  
>   
>     
>       
> X
>       
>     
> {\\displaystyle X}
>   
> and 
>   
>     
>       
> Y
>       
>     
> {\\displaystyle Y}
>   
> , are the square roots of the (necessarily non-negative) eigenvalues of the self-adjoint operator 
>   
>     
>       
>         
> ...

“ The singular values are the absolute values of the [eigenvalues](https://en.wikipedia.org/api/rest_v1/page/mobile-html/Eigenvalues) of a [normal matrix](https://en.wikipedia.org/api/rest_v1/page/mobile-html/Normal_matrix) _A_ , because the [spectral theorem](https://en.wikipedia.org/api/rest_v1/page/mobile-html/Spectral_theorem) can be applied to obtain unitary diagonalization of ![A](https://wikimedia.org/api/rest_v1/media/math/render/svg/7daff47fa58cdfd29dc333def748ff5fa4c923e3) as ![{isplaystyle A=Uambda U^{*}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3f8ccf13ba3314ff6981a78452b2e23b76998d5d). Therefore, ![{extstyle {qrt {A^{}A}}={qrt {Uambda ^{}ambda U^{}}}=Ueft|ambda ight|U^{}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e0db0f581b7bc0239ec3914febece3d7aa3099d2).”

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [May 4, 2023, 10:36am UTC](https://discourse.julialang.org/t/absolute-value-of-a-matrix/98301/3 "2023-05-04T10:36:03Z")

</div>

The first two methods would not stable in general if the matrix has singular values of widely varying magnitudes. You narrowed the focus to Hermitian matrices, but this is defined more generally and `absmat3` is numerically fine for Hermitian matrices, but it is wrong for general square matrices.

More generally, \sqrt{A^\* A} is the semidefinite factor H in the polar factorization A=UH where U is unitary and can be obtained from the SVD:

```julia
function absmat4(A::AbstractMatrix)
    F = svd(A)
    return F.Vt' * Diagonal(F.S) * F.Vt
end

```

That is not faster, but works for any square matrix. You could then have

```julia
function absmat3(A::Hermitian{<:AbstractMatrix})
    F = eigen(A)
    return F.vectors * Diagonal(abs.(F.values)) * F.vectors'
end

```

I don’t think you will do much better in terms of methods based on the standard matrix decompositions. But there are iterative methods for U which could then give you H=U^\* A. A Newton iteration for computing U takes the form

X\_{k+1} = \frac{1}{2} (X\_k + X\_k^{-\*}), \qquad X\_0 = A.

X\_k is quadratically convergent to U. There also appears to be a package for this already: [PolarFact.jl](https://github.com/weijianzhang/PolarFact.jl). The README provides some references.

---

<div class="post-metadata">

**Author:** ![platawiec](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/platawiec/32/31914_2.png) [@platawiec](https://discourse.julialang.org/u/platawiec)\
**Post date:** [May 4, 2023, 11:59am UTC](https://discourse.julialang.org/t/absolute-value-of-a-matrix/98301/4 "2023-05-04T11:59:27Z")

</div>

Great answer. A (non iterative) polar factorization is also available here [GitHub - JuliaLinearAlgebra/MatrixFactorizations.jl: A Julia package to contain non-standard matrix factorizations](https://github.com/JuliaLinearAlgebra/MatrixFactorizations.jl)