# Inverse of symbolic matrix

``````using Symbolics, LinearAlgebra
@variables q::Real
iden = collect(I(4) .* 1.0)
k = inv(iden .* q^2)
``````

I get true/q^2 on the diagonal instead of 1/q^2

2 Likes
``````julia> true == 1
true
``````

`true` is the smallest value that represents the integer `1`

``````julia> true isa Integer
true
``````
1 Like

This looks like a bug that is worth reporting. Even more extreme is

``````using Symbolics, LinearAlgebra
@variables q
inv(collect(q^2*2.0*I(4)))
``````

giving

``````4×4 Matrix{Num}:
true / (2.0(q^2))    0                  0                  0
0                true / (2.0(q^2))    0                  0
0                  0                true / (2.0(q^2))    0
0                  0                  0                true / (2.0(q^2))
``````
1 Like

What do you find problematic? What answer would you have expected, taking into account that `true == 1`?

1 Like

With more complicated expressions it becomes difficult to read 1 instead of true. I tried “simplify” but it doesn’t convert true to 1. Maybe there is an alternative? I am declaring q as Real and would not want a bool result.

There is only a `Bool` in the symbolic expression, the result will be a floating-point value.

``````julia> q = rand();

julia> true/2q^2 === 1/2q^2 # Absolutely identical
true

julia> true/2q^2 # Result is a floating-point value
2.8977550256308775

julia> true/2q^2 |> typeof
Float64
``````

The whole point is that `true == 1`, just like `1.0 == 1` , and it makes sense to use the smallest version of 1 possible to avoid type promotion if not required.

1 Like

You’re right, but I might want to use Symbolics to get a symbolic result without necessarily wanting to calculate it numerically.
In a symbolic expression I would like numbers, not equivalent quantities like Bool. Thank you.

1 Like

You could modify how `Bool` is printed, but it may have unintended consequences

``````julia> Base.show(io::IO, b::Bool) = print(io, (b ? "1" : "0"))

julia> true
1

julia> false
0
``````
1 Like

But they are not equivalent. `1.0 === 1` is false and we can get different answers in subsequent computations

``````julia> x = sum(1.0 for i in range(1,2^16))
65536.0

julia> 2^48*x
1.8446744073709552e19

julia> x = sum(1 for i in range(1,2^16))
65536

julia> 2^48*x
0
``````
1 Like

If you write the prompt more concisely, the result is more readable, too.

``````julia> using Symbolics, LinearAlgebra

julia> @variables q
1-element Vector{Num}:
q

julia> inv(collect(q^2*2I(4)))
4×4 Matrix{Num}:
1 / (2(q^2))             0             0             0
0  1 / (2(q^2))             0             0
0             0  1 / (2(q^2))             0
0             0             0  1 / (2(q^2))

``````
1 Like

The computations you show here are unrelated to the expression `true/2q^2` which is in fact identical to `1/2q^2`. If you intend to make further numerical computations, you can always make sure you have a suitable type for the computations you want to perform.

1 Like

To elaborate on why `true` is a good choice here, consider the following trivial function, if we define it using `false`, we see that `false` is adaptive in the sense that multiplying with `false` will make the result have the same type as the input, irrespective of what the type is:

``````julia> multiply_by_zero(x) = false*x
multiply_by_zero (generic function with 1 method)

julia> multiply_by_zero(true)
false

julia> multiply_by_zero(1)
0

julia> multiply_by_zero(1f0)
0.0f0

julia> multiply_by_zero(1.0)
0.0

julia> multiply_by_zero(big"1.0")
0.0
``````

Had we defined the function using a floating-point value, we see that the type of the result is different from the type of the input whenever the input type is smaller than `Float64`

``````julia> multiply_by_zero(x) = 0.0*x
multiply_by_zero (generic function with 1 method)

julia> multiply_by_zero(true)
0.0

julia> multiply_by_zero(1)
0.0

julia> multiply_by_zero(1f0)
0.0

julia> multiply_by_zero(1.0)
0.0

julia> multiply_by_zero(big"1.0")
0.0

julia>
``````
3 Likes

With a float coefficient that is unfortunately not the case

``````julia> using Symbolics, LinearAlgebra

julia> @variables q
1-element Vector{Num}:
q

julia> z = inv(collect(q^2*2.0I(4)))
4×4 Matrix{Num}:
true / (2.0(q^2))    0                  0                  0
0                true / (2.0(q^2))    0                  0
0                  0                true / (2.0(q^2))    0
0                  0                  0                true / (2.0(q^2))
``````

Going further with this example worries me

``````julia> y = z[1,1]
true / (2.0(q^2))

julia> x = sum(y for i in range(1,2^16))
65536(true / (2.0(q^2)))

julia> 2^48*x
0(true / (2.0(q^2)))
``````
1 Like

simplifying or expanding doesn’t seem to help either:

``````julia> simplify(y)
true / (2.0(q^2))

julia> simplify(x)
65536(true / (2.0(q^2)))

julia> expand(x)
65536 / (2.0(q^2))

julia> 2^48*expand(x)
0

julia> simplify(expand(x))
65536 / (2.0(q^2))
``````
1 Like

This is simply integer overflow

``````julia> 2^48
281474976710656

julia> 65536 * 2^48
0
``````

it is completely unrelated to Symbolics

1 Like

qed

``````julia> 2^48
281474976710656

julia> 65536 * 2^48
0

julia> 65536 * BigInt(2)^48
18446744073709551616

``````
1 Like

Yes, that is how I constructed the example. If you reduce your float to an integer then you run in to the overflow problem in subsequent operations. If you keep it is a float, then that problem disappears, at the expense of the answer not being exactly representable. If you change 2.0 to 2.1 or 2.01063072 in the above examples the problems persist, and in the latter cases you would not expect simple replacements with integral or rational constants.

How does one get symbolics to reduce

``````true/(2.0*q^2)
``````

to a more consistent form such as

``````0.5 q^(-2)
``````

or

``````0.5 / q^2
``````

?

2 Likes

I see what you mean; the overflow is hard to diagnose if it accumulates in different parts of the code.

First option I can think of: Use `BigInt` everywhere; that should be okay in a symbolic computation.

Second option is to implement a `float` method which recursively turns all numeric coeffiecients into floating-point numbers. see Symbolics.jl/num.jl at master · JuliaSymbolics/Symbolics.jl · GitHub

1 Like

The Num type has builtin choices for its zero element and its unit element given by

``````julia> using Symbolics

julia> zero(Num)
0

julia> typeof(zero(Num))
Num

julia> one(Num)
1

julia> typeof(one(Num))
Num
``````

We have for example

``````julia> one(Num)/2.0
0.5

julia> one(Num)/(2.0*q^2)
1 / (2.0(q^2))
``````

So why is this not used?

2 Likes

I would have expected either

``````0.5 / q^2
``````

or

``````0.5 q^(-2)
``````

rather than

``````true / (2.0(q^2))
``````
1 Like