# Cholesky promotes Float16

**URL:** <https://discourse.julialang.org/t/cholesky-promotes-float16/76772>\
**Category:** Numerics\
**Tags:** question\
**Created:** [February 19, 2022, 7:07pm UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772 "2022-02-19T19:07:07Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [February 19, 2022, 7:07pm UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/1 "2022-02-19T19:07:08Z")

</div>

I’m seeing Cholesky promote half precision to single on its own. Is this  
supposed to happen? lu does not do this.

Running 1.7.2 on an M1 Macbook. Same results on Intel Mac with MKL Blas.

Bug or feature?

```julia
julia> R=rand(Float16,3,3)
3×3 Matrix{Float16}:
 4.00391e-01 8.24219e-01 2.50000e-01
 2.34375e-02 2.90527e-01 1.09375e-01
 7.85156e-01 9.42383e-01 3.76465e-01

julia> P=R'*R
3×3 Matrix{Float16}:
 7.77344e-01 1.07617e+00 3.98438e-01
 1.07617e+00 1.65234e+00 5.92773e-01
 3.98438e-01 5.92773e-01 2.16187e-01

julia> PF=lu(P)
LU{Float16, Matrix{Float16}}
L factor:
3×3 Matrix{Float16}:
 1.00000e+00 0.00000e+00 0.00000e+00
 7.22168e-01 1.00000e+00 0.00000e+00
 3.70117e-01 1.58325e-01 1.00000e+00
U factor:
3×3 Matrix{Float16}:
 1.07617e+00 1.65234e+00 5.92773e-01
 0.00000e+00 -1.17188e-01 -2.95410e-02
 0.00000e+00 0.00000e+00 1.50299e-03

julia> QF=cholesky(P)
Cholesky{Float32, Matrix{Float32}}
U factor:
3×3 UpperTriangular{Float32, Matrix{Float32}}:
 8.81671e-01 1.22060e+00 4.51912e-01
      ⋅ 4.03073e-01 1.02135e-01
      ⋅ ⋅ 3.91250e-02

```

---

<div class="post-metadata">

**Author:** ![artkuo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/artkuo/32/28577_2.png) [@artkuo](https://discourse.julialang.org/u/artkuo)\
**Post date:** [February 19, 2022, 8:02pm UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/2 "2022-02-19T20:02:32Z")

</div>

It looks like it’s done on purpose. It promotes to at least `Float32`, based on the type of the square root of an element of the matrix:  
[https://github.com/JuliaLang/julia/blob/623ceb7834de47538eddeadfa84a8bf2d9741248/stdlib/LinearAlgebra/src/cholesky.jl#L181](https://github.com/JuliaLang/julia/blob/623ceb7834de47538eddeadfa84a8bf2d9741248/stdlib/LinearAlgebra/src/cholesky.jl#L181)  
I’m not sure why, but it could be that they wanted to ensure some minimum precision for the result. It might be worth debating whether that promotion should be made optional, or at least document that it does so.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [February 19, 2022, 8:16pm UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/3 "2022-02-19T20:16:08Z")

</div>

Nice work to find that.

My view is that lu does the right thing. If you want to compute in Float16, the language should let you. I think this is likely an unintended consequence of that line. It would  
seem that

`@inline choltype(A) = typeof(sqrt(oneunit(eltype(A)))`  
would make more sense.

However, I’m not clear that it would not break something else.

---

<div class="post-metadata">

**Author:** ![artkuo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/artkuo/32/28577_2.png) [@artkuo](https://discourse.julialang.org/u/artkuo)\
**Post date:** [February 19, 2022, 9:17pm UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/4 "2022-02-19T21:17:54Z")

</div>

I can make a guess for the reason, which makes this a bug. In the `cholesky.jl` code,

```julia
# The dispatch structure in the cholesky, and cholesky! methods is a bit
# complicated and some explanation is therefore provided in the following
#
# In the methods below, LAPACK is called when possible, i.e. StridedMatrices with Float32,
# Float64, ComplexF32, and ComplexF64 element types. 

```

My guess it that `choltype` was originally intended to ensure at least Float32 for the LAPACK branch, but is accidentally used even for the types that don’t get sent to LAPACK. In fact, later in the same comments,

```julia
# FixMe? The dispatch below seems overly complicated. One simplification could be to
# merge the two Cholesky types into one. 

```

So I believe the correct fix would be to keep both yours and the original version (e.g. by naming `choltypeforlapack` or by dispatch). As you say, it’s not clear if your fix would break something else, but I would label the current non-LAPACK behavior a bug.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [February 19, 2022, 9:51pm UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/5 "2022-02-19T21:51:03Z")

</div>

I raised an issue. We’ll see what the complier people think.

---

<div class="post-metadata">

**Author:** ![gbaraldi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gbaraldi/32/22101_2.png) [@gbaraldi](https://discourse.julialang.org/u/gbaraldi)\
**Post date:** [February 20, 2022, 1:47am UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/6 "2022-02-20T01:47:57Z")

</div>

I think there are two issues at play. First thing is that the LAPACK implementation that julia uses doesn’t have float16 routines which is probably why it has to get promoted, the other is that julia currently has a compiler pass to promote Float16 computations to Float32.  
One thing that could be done is to truncate the result after doing the LAPACK routine so the user gets a Float16 back.  
I think in a previous thread you opened I gave more details [Half precision](https://discourse.julialang.org/t/half-precision/75580)

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [February 20, 2022, 12:42pm UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/7 "2022-02-20T12:42:09Z")

</div>

Strange. Everything works as I would expect with `lu`. In particular, `lu!` takes a long time in Float16 because Float16 computations are done in software.

```julia
julia> A=rand(Float32, 1000,1000); A16=Float16.(A);

julia> @btime lu!($A);
  4.907 ms (1 allocation: 7.94 KiB)

julia> @btime lu!($A16);
  218.729 ms (1 allocation: 7.94 KiB)

```

However, `cholesky` seems to do all of its work in Float32, so the matrix gets promoted on the way in before LAPACK. You can see that the timings are much closer than for `lu`. That promotion makes some sense given LAPACK’s limits. Why is `lu` different? Somehow I can’t get `cholesky!` to work right with `@btime`.

```julia
julia> A=rand(Float32,1000,1000); A ./= opnorm(A);

julia> B=I + Hermitian(A*A); B16=Float16.(B);

julia> @btime cholesky($B);
  3.025 ms (4 allocations: 3.81 MiB)

julia> @btime cholesky($B16);
  3.386 ms (4 allocations: 3.81 MiB)

```

---

<div class="post-metadata">

**Author:** ![gbaraldi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gbaraldi/32/22101_2.png) [@gbaraldi](https://discourse.julialang.org/u/gbaraldi)\
**Post date:** [February 20, 2022, 12:56pm UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/8 "2022-02-20T12:56:44Z")

</div>

Probably LU has a simple fallback implementation in Julia while cholesky doesn’t.

---

<div class="post-metadata">

**Author:** ![dkarrasch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dkarrasch/32/7410_2.png) [@dkarrasch](https://discourse.julialang.org/u/dkarrasch)\
**Post date:** [February 20, 2022, 5:35pm UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/9 "2022-02-20T17:35:08Z")

</div>

The type promotion to at least Float32 is indeed to make use of LAPACK. Converting back down to Float16 was implemented in [https://github.com/JuliaLang/julia/pull/41352](https://github.com/JuliaLang/julia/pull/41352) and is going to be available in v1.8.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [February 20, 2022, 6:15pm UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/10 "2022-02-20T18:15:02Z")

</div>

Just what I wanted to hear. Thanks.

---

<div class="post-metadata">

**Author:** ![Paul\_Soderlind](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paul_soderlind/32/1753_2.png) [@Paul\_Soderlind](https://discourse.julialang.org/u/Paul_Soderlind)\
**Post date:** [February 20, 2022, 7:30pm UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/11 "2022-02-20T19:30:56Z")

</div>

A related issue is that `lowrankupdate!(C,v)` fails if `v` is a vector of integers.

---

<div class="post-metadata">

**Author:** ![artkuo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/artkuo/32/28577_2.png) [@artkuo](https://discourse.julialang.org/u/artkuo)\
**Post date:** [February 20, 2022, 7:31pm UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/12 "2022-02-20T19:31:22Z")

</div>

That’s good news. But @ctkelley note that this still promotes and computes in `Float32` and then converts the output down to `Float16`. It still promotes in `choltype`, whether for LAPACK or generic. In this sense it’s still different from `lu`. It does ensures the right type is returned, but if what you want is a true `Float16` emulation, the precision here may actually be a bit higher.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [February 20, 2022, 8:10pm UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/13 "2022-02-20T20:10:32Z")

</div>

The main reason for using `Float16` and such is performance.  
Converting it to `Float32`, paying the run time tax and then going back is the worse.  
Something like not having the cake and not eating it.

If it is an intermediate phase as the current implementations (_LAPACK_) doesn’t support `Float16` it is OK. But the end goal should be that using `Float16` must be faster than `Float32` and _might be_ less accurate.

Is there a plan for native Linear Algebra wide support for `Float16`?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [February 20, 2022, 8:13pm UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/14 "2022-02-20T20:13:35Z")

</div>

float16 can’t be faster for most processors since most don’t have hardware float16 support.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [February 20, 2022, 9:27pm UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/15 "2022-02-20T21:27:20Z")

</div>

Apple silicon supports Float16, as do the leadership class supercomputers. There is some pretty wide and serious interest in this.

> **[A Survey of Numerical Methods Utilizing Mixed Precision Arithmetic](https://arxiv.org/abs/2007.06674)**
>
> Within the past years, hardware vendors have started designing low precision special function units in response to the demand of the Machine Learning community and their demand for high compute power in low precision formats. Also the server-line...

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [February 21, 2022, 2:40am UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/16 "2022-02-21T02:40:20Z")

</div>

See [https://github.com/JuliaLang/julia/issues/40216#issuecomment-1046411397](https://github.com/JuliaLang/julia/issues/40216#issuecomment-1046411397) for progress on getting better FP16 support for Julia.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [February 21, 2022, 5:24am UTC](https://discourse.julialang.org/t/cholesky-promotes-float16/76772/17 "2022-02-21T05:24:53Z")

</div>

You forget one important element for performance - Throughput.  
While the execution ports of most CPU’s doesn’t support `Float16` natively, for most Linear Algebra operation the issue isn’t the rate of execution but the memory bandwidth which is effectively multiplied (Unless in memory in Julia array of `Float16` isn’t dense).

But what I wrote was the end goal. We’re also during intermediate phase in CPU’s as well.  
But what I wrote would have happened on CPU’s with `Float16` support, so it holds.