JuliaReach
Symbolic interval matrix power type
This issue is related to https://github.com/JuliaReach/IntervalMatrices.jl/issues/88. Proposal: add a new IntervalMatrixPower wrapper which holds the power k of the matrix power. The purpose of this type is to experiment with lazy construction of the matrix power for evaluation using different methods. The new type can also have a cache for the representation of the power using symbols, eg. polynomials.
I agree. In the application we only want to increment the power symbolically, i.e., from A^k we want to compute A^{k+1}. Proposal (not tested):
struct IntervalMatrixPower
M::SymbolicMatrix # `SymbolicMatrix` is the type of a symbolic matrix
Mᵏ::SymbolicMatrix
function IntervalMatrixPower(M::IntervalMatrix, k::Int)
Ms = symbolic(M) # `symbolic` converts to a symbolic matrix
Msᵏ = Ms^k
return new(Ms, Msᵏ)
end
end
increment!(Aᵏ::IntervalMatrixPower) = Aᵏ.Mᵏ *= Aᵏ.M
increment(Aᵏ::IntervalMatrixPower) = increment!(copy(Aᵏ))
Possible problem: k is not stored explicitly. That would require a mutable struct.
True, that implementation is more or less what i had in mind. It may be convenient to subtype <: AbstractIntervalMatrix. But this requires more work on accessing the coefficients, etc.
I think that the true gain will come from a tight evaluation of the polynomial coefficients of M^k. If we just use the result of M^k computed symbolically that is not going to help too much. Using the Horner form + interval arithmetic would help. Otherwise, we can use tighter enclosure methods with RangeEnclosures.
Perhaps before investing too much on this, we can consider a reachability example manually and see if there are considerable gains with using M^k in other way that the naive matrix multiplication.
It may be convenient to subtype <: AbstractIntervalMatrix. But this requires more work on accessing the coefficients, etc.
I would keep this for later. Such an implementation is definitely more time consuming and I do not see the benefit (we need specific code in the application anyway).
If we just use the result of M^k computed symbolically that is not going to help too much.
I believe that knowing the dependencies can already gain a lot, assuming that the concrete interval power is clever (i.e., if x^k is more precise than k multiplications of x). I do not understand how the "tight evaluation" would look like.
Perhaps before investing too much on this, we can consider a reachability example manually and see if there are considerable gains with using M^k in other way that the naive matrix multiplication.
Why reachability? Just evaluate A^k explicitly and symbolically.
Why reachability? Just evaluate A^k explicitly and symbolically.
True, that shows a difference, but i don't know how much of that difference we need to observe a meaningful consequence in reachability (im thinking in the transmission line model). It depends on the actual numerical coefficients of the model.
Here is code to compare between the evaluation of the symbolic matrix (without post-processing) vs. normal interval matrix multiplication (the full code is in this notebook):
# self-contained example
using SymEngine, RangeEnclosures, IntervalMatrices
using MacroTools: postwalk
subidx(i) = join(Char.(0x2080 .+ convert.(UInt16, digits(i)[end:-1:1])))
function symbolic_matrix(n; name="M")
M = Matrix{SymEngine.Basic}(undef, n, n)
for i in 1:n, j in 1:n
M[i, j] = name * subidx(i) * subidx(j)
end
return M
end
function subs(ex, imat; name="M")
n = size(imat, 1) # imat assumed square
for i in 1:n, j in 1:n
name_ij = Symbol(name * subidx(i) * subidx(j))
ex == name_ij && return imat[i, j]
end
return ex
end
julia> M1 = rand(IntervalMatrix) # M^1 as an interval matrix (known exactly)
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
[-0.988735, -0.0794308] [-1.31914, -1.16598]
[-0.467652, 2.06963] [-0.90236, 0.14924]
julia> M = symbolic_matrix(2)
2×2 Array{Basic,2}:
M₁₁ M₁₂
M₂₁ M₂₂
julia> M2 = M * M
2×2 Array{Basic,2}:
M₂₁*M₁₂ + M₁₁^2 M₁₁*M₁₂ + M₂₂*M₁₂
M₂₁*M₁₁ + M₂₂*M₂₁ M₂₁*M₁₂ + M₂₂^2
julia> Msym(k) = M^k
julia> Msym_eval(k) = [postwalk(ex -> subs(ex, M1), convert(Expr, m)) for m in Msym(k)];
julia> eval.(Msym_eval(10))
2×2 Array{Interval{Float64},2}:
[-1738.29, 1817.95] [-1448.11, 1354.38]
[-2127.52, 2274.61] [-1739.91, 1819.78]
julia> M1^10
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
[-2517.86, 2293.8] [-2182.93, 1652.21]
[-2680.61, 3477.52] [-2549.71, 2380.51]
I think the wrapper should combine #83 and #88: Every square number can and should be computed exactly (and rather efficiently). (Note that this is more precise than the symbolic evaluation.) However, there are not too many of them :cry:
For the numbers in between, the symbolic version can be used. Unfortunately, the symbolic version cannot make use of the intermediate results (i.e., the square powers) (that was the idea in #83). But I think we want to support different ways of evaluating the power because we saw that the precise evaluation is expensive. For example, there can be eval!(p, algorithm="symbolic") for the precise evaluation and eval!(p, algorithm="fast") which just takes the previous result and multiplies it with M.
Continuing my example from above, IntervalMatrixPower would also store M_concrete and Mᵏ_concrete and k and k_cached):
function eval!(p::IntervalMatrixPower; algorithm::String="fast")
if p.k_cached == p.k
return p.Mᵏ_concrete
elseif isapoweroftwo(p.k)
result = eval_poweroftwo(p)
elseif algorithm == "symbolic"
result = eval_symbolic(p)
elseif algorithm == "fast"
result = eval_fast(p)
else
throw(ArgumentError("algorithm $algorithm unknown"))
end
p.Mᵏ_concrete = result
p.k_cached = p.k
return result
end
function isapoweroftwo(k::Integer)
(k & (k - 1)) == 0 # see https://stackoverflow.com/a/600306
end
function eval_poweroftwo(p::IntervalMatrixPower)
Mⁱ = copy(p.M_concrete)
@inbounds for i in 1:Int(log(2, p.k))
Mⁱ = square(Mⁱ)
end
return Mⁱ
end
function eval_symbolic(p::IntervalMatrixPower)
... # as before
end
function eval_fast(p::IntervalMatrixPower)
return p.Mᵏ_concrete * p.M_concrete
end