# Wrong output from det(Array{Complex})

**URL:** https://discourse.julialang.org/t/wrong-output-from-det-array-complex/10022
**Category:** General Usage
**Tags:** linearalgebra
**Created:** [March 27, 2018, 10:18pm UTC](https://discourse.julialang.org/t/wrong-output-from-det-array-complex/10022 "2018-03-27T22:18:20Z")
**Posts on this page:** 12
**Page:** 1

<div class="post-metadata">

### Author: ![pedrohnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pedrohnv/32/212503_2.png) [@pedrohnv](https://discourse.julialang.org/u/pedrohnv)
#### Post date: [March 27, 2018, 10:18pm UTC](https://discourse.julialang.org/t/wrong-output-from-det-array-complex/10022/1 "2018-03-27T22:18:20Z")

</div>

I have a complex matrix

`A = Array{Complex32}([20-50im -10+20im -10+30im; -10.0+20.0im 26.0-52.0im -16.0+32.0im; -10.0+30.0im -16.0+32.0im 26.0-62.0im])`

Its determinant is 0, but when I call

```julia
det(A) #out=-116.25 - 5.00im, should be 0
det(real(A)) #out=0.0, correct
det(imag(A)) #put=-68.75, should also be 0

```

Strangely, when I promote the imaginary part, it gives the correct result

`det(Array{Real}(imag(A))) # outputs 0.0`

I suppose it has to do with the precision used to create A (Complex32). When I create it using Complex128, the determinant of the real and imaginary part are correct, but not of the matrix as a whole.

```julia
A = Array{Complex64}([20-50im -10+20im -10+30im; -10.0+20.0im 26.0-52.0im -16.0+32.0im; -10.0+30.0im -16.0+32.0im 26.0-62.0im])
det(A) #-0.0074005127 + 0.006790161im
det(real(A)) #0.0
det(imag(A)) #0.0

```

My question is, why is it that det(A) gives the wrong result (error too big, that is) and how can I impove it?

---

<div class="post-metadata">

### Author: ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)
#### Post date: [March 27, 2018, 11:33pm UTC](https://discourse.julialang.org/t/wrong-output-from-det-array-complex/10022/2 "2018-03-27T23:33:26Z")

</div>

Note that `Complex32` stores its real and imaginary components as `Float16`s (so that the total size of a `Complex32` is 32 bytes). `Float16` is a _really_ small floating point type, so large errors can be expected. For `Complex128`, the determinant is `-4.654054919228656e-11 + 4.618527782440645e-12im`, which isn’t too bad. Floating point math is just different from math with real numbers, so you may have to lower your expectations a bit when using them. See also this excellent post: [PSA: floating-point arithmetic - #5 by Tamas\_Papp](https://discourse.julialang.org/t/psa-floating-point-arithmetic/8678/5).

If you really need a lot of precision (at the cost of performance), you could use `Complex{BigFloat}`. If you want to go that route, note that something like `26.0-52.0im` constructs a `Complex{Float64}`, and if the number you want isn’t exactly representable as a `Complex{Float64}`, you’ll lose precision already before converting to `BigFloat`. Alternatively, in your specific case, the real and imaginary parts are integers, so if you replace e.g. `26.0-52.0im` with `26-52im` and construct an `Array{BigFloat}`, you won’t lose precision.

See also: [Accuracy of Computing Determinants - #3 by andreasnoack](https://discourse.julialang.org/t/accuracy-of-computing-determinants/6782/3).

---

<div class="post-metadata">

### Author: ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 27, 2018, 11:59pm UTC](https://discourse.julialang.org/t/wrong-output-from-det-array-complex/10022/3 "2018-03-27T23:59:32Z")

</div>

> [@pedrohnv](#):
>
> My question is, why is it that `det(A)` gives the wrong result (error too big, that is) and how can I impove it?

First of all, you should understand that accurate computation of the determinant is a notoriously difficult problem. See, for example,

- Erich Kaltofen and Gilles Villard, [“Computing the sign or the value of the determinant of an integer matrix, a complexity survey,”](https://www.sciencedirect.com/science/article/pii/S0377042703007234) _J. Comp. Appl. Math_ 162, pp. 133–146 (2004).

In particular, that paper says that the condition number of `det(A)` for an n×n matrix `A` is on the order of `norm(A, Inf)^n`. For your matrix, that gives about 2.44×10⁶. So, I would naively expect that computing determinants of your matrix could easily lose 6 digits to roundoff errors, which means that for `Float16` arithmetic (which keeps 3 digits) I would expect to lose everything, and even `Float32` arithmetic (which keeps about 7 digits) is problematic.

`det(A)` works by performing the LU factorization of `A` and then multiplying the diagonal entries of `U`. This is efficient for large n, O(n³), but it incurs roundoff errors even for integer-valued matrix entries.

For a 3×3 matrix, you can instead use the direct permutation formula for the determinant, which is exact for integer values that aren’t too big. This is built-in to the StaticArrays.jl package (it is a good idea to use StaticArrays anyway, for efficiency, if you are working with fixed size ≤ 4×4 matrices):

```julia
julia> using StaticArrays

julia> A = @SMatrix Complex{Float32}[20-50im -10+20im -10+30im; -10.0+20.0im 26.0-52.0im -16.0+32.0im; -10.0+30.0im -16.0+32.0im 26.0-62.0im];

julia> det(A)
0.0f0 + 0.0f0im

```

Even with `SMatrix`, `Float16` won’t work because the individual terms in the determinant formula overflow the maximum `Float16` value. But it is **really inadvisable** to do computations with `Float16` anyway — `Float32` or `Float64` are much faster. `Float16` is mainly useful for low-precision _storage_ of the result.

All this being said, you should really **think carefully about whether you need to compute the determinant at all.** It is rarely used nowadays in numerical computations — there is often a better way to do whatever you wanted the determinant for.

---

<div class="post-metadata">

### Author: ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [March 28, 2018, 6:46am UTC](https://discourse.julialang.org/t/wrong-output-from-det-array-complex/10022/4 "2018-03-28T06:46:15Z")

</div>

In addition to the excellent suggestions of @tkoolen and @stevengj, I would look at the context: _why_ do you need determinants? If you are exploring the stucture of a matrix (eg rank), QR, or ideally SVD are usually better options for numerical linear algebra. Using 64 bit floats (ie `Complex{Float64}` in your example) are usually worth it.

---

<div class="post-metadata">

### Author: ![pedrohnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pedrohnv/32/212503_2.png) [@pedrohnv](https://discourse.julialang.org/u/pedrohnv)
#### Post date: [March 28, 2018, 1:42pm UTC](https://discourse.julialang.org/t/wrong-output-from-det-array-complex/10022/5 "2018-03-28T13:42:16Z")

</div>

Thank you for your reply,

Actually, I do not need the determinant for my problem. I was just toying with julia. But I do need to solve a big linear algebra problem, though I just need approximate results (any error \< 1e-2 is acceptable).

---

<div class="post-metadata">

### Author: ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [March 28, 2018, 1:47pm UTC](https://discourse.julialang.org/t/wrong-output-from-det-array-complex/10022/6 "2018-03-28T13:47:19Z")

</div>

Depending on the problem, iterative methods may be a good choice for large problems in need of a “good enough” solution. Try  
[https://github.com/JuliaMath/IterativeSolvers.jl](https://github.com/JuliaMath/IterativeSolvers.jl)

---

<div class="post-metadata">

### Author: ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 28, 2018, 3:54pm UTC](https://discourse.julialang.org/t/wrong-output-from-det-array-complex/10022/7 "2018-03-28T15:54:49Z")

</div>

If you want to solve Ax=b, look at the accuracy of `A \ b` (which depends on the conditioning of A), not determinants.

Whether to do something fancier like iterative solvers depends on what you mean by “big”. (1000x1000 is not big, for example.)

---

<div class="post-metadata">

### Author: ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [March 30, 2018, 3:20pm UTC](https://discourse.julialang.org/t/wrong-output-from-det-array-complex/10022/8 "2018-03-30T15:20:10Z")

</div>

> [@tkoolen](#):
>
> Note that Complex32 stores its real and imaginary components as Float16s (so that the total size of a Complex32 is 32 bytes). Fl

Are you sure? I think 32 stands for the number of bits in both the real and imaginary parts.

---

<div class="post-metadata">

### Author: ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)
#### Post date: [March 30, 2018, 3:27pm UTC](https://discourse.julialang.org/t/wrong-output-from-det-array-complex/10022/9 "2018-03-30T15:27:34Z")

</div>

```julia
julia> dump(Complex32)
Complex{Float16} <: Number
  re::Float16
  im::Float16

```

---

<div class="post-metadata">

### Author: ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [March 30, 2018, 4:05pm UTC](https://discourse.julialang.org/t/wrong-output-from-det-array-complex/10022/10 "2018-03-30T16:05:50Z")

</div>

Right, one letter difference: I confused it with `ComplexF32`…

---

<div class="post-metadata">

### Author: ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)
#### Post date: [March 30, 2018, 5:15pm UTC](https://discourse.julialang.org/t/wrong-output-from-det-array-complex/10022/11 "2018-03-30T17:15:13Z")

</div>

Ah, I guess this will be changed in 0.7 (`ComplexF32` is not a thing in 0.6):

```julia
julia> ComplexF32
Complex{Float32}

julia> Complex32
WARNING: Base.Complex32 is deprecated, use ComplexF16 instead.
 in module Main
Complex{Float16}

```

---

<div class="post-metadata">

### Author: ![greg\_plowman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/greg_plowman/32/8100_2.png) [@greg\_plowman](https://discourse.julialang.org/u/greg_plowman)
#### Post date: [March 31, 2018, 2:54am UTC](https://discourse.julialang.org/t/wrong-output-from-det-array-complex/10022/12 "2018-03-31T02:54:14Z")

</div>

> [@PetrKryslUCSD](#):
>
> Are you sure? I think 32 stands for the number of bits in both the real and imaginary parts.

> [@tkoolen](#):
>
> Ah, I guess this will be changed in 0.7 (ComplexF32 is not a thing in 0.6)

I think this has long been a source of ambiguity and confusion.

See discussion and PR:

> [@Rename Complex2n to Complexn?](https://discourse.julialang.org/t/rename-complex2n-to-complexn/6081):
>
> I have always been bugged by the names Complex128 = Complex{Float64}, etc, as there is always (probably because I’m not a heavy user) a tiny mental step to divide the number to know the actual precision (64 bits for Complex128). I would find it clearer to have Complex64 meaning Complex{Float64}, so that the number indicates directly the float type used. Also, it seems more scalable: Quaternion256 starts to be unreadable, let alone Octonion512. My concern is also that it makes a precedent for nam…

> <https://github.com/JuliaLang/julia/pull/24647>
>
> Cf. https://discourse.julialang.org/t/rename-complex2n-to-complexn/6081 for back…ground.
> Here is a recollection from this thread of 3 reasons to deprecate the \`ComplexN\` aliases:
> 
> \- In \`Complex128\`, \`128\` refers to no semantic reality connected to the type. It only refers to the very low-level fact that it's its number of bits. Actually, in the initial discourse thread, I suggested renaming \`Complex128\` to \`Complex64\` with the idea "the coordinates of a \`Complex64\` is \`Float64\`", so there is a clearer connection between the 2 types.
> \- \`Complex128\` is ambiguous, since any eltype with 64 bit size (e.g. Int64) will give you a Complex type of size 128 (@Keno);
> \- Historically, the \`Complex128\` type actually preceded \`Complex{Float64}\` — it was originally a bits type because Julia didn’t have immutable types. So, at this point, deprecating \`Complex128\` makes a lot of sense (@stevengj).
