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