Symbolic interval matrix power type

[mforets]

This issue is related to #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.

[schillic]

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.

[mforets]

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.

[schillic]

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.

[mforets]

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.

[mforets]

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]

[schillic]

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 😢

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