Matrix power via square decomposition

[schillic]

Problem

The k-th integer power of an interval matrix M, M^k, is not well-defined. The reason is that interval-matrix multiplication is not associative. Furthermore, even for a fixed order, the multiplication induces an error due to the intervals. The only safe case is squaring (i.e., k = 2) (see #79).

Idea

Decompose the power k such that many squares are used. In other words: use as few non-square multiplications as possible.

Example

In the example below, we compute M^9.

• A and B use the naive power operation, which are the worst options. (Interestingly, ^ does something slightly more intelligent. ?> ^ says equivalent to \exp(p\log(A))).
• C, D, E use the fact that 9 = 2² + 2² + 1.
• F, G, H use the fact that 9 = (2 + 1)³ + 2 + 1 (I only tried three permutations).
• As can be seen, the order of the multiplications plays a role.

julia> M = rand(IntervalMatrix)
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-1.88661, 1.81318]     [-0.483837, 0.892705]
 [-0.893564, -0.447015]  [-0.10882, 1.85449]

julia> A = M*M*M*M*M*M*M*M*M  # 9 multiplications
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-4394.68, 4388.85]  [-4300.39, 4320.52]
 [-3599.78, 3582.38]  [-3495.86, 3545.53]

julia> B = M^9
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-3734.67, 3685.31]  [-3311.99, 3430.62]
 [-3602.63, 3458.46]  [-3008.63, 3290.32]

julia> M2 = square(M); M4 = square(M2);

julia> C = M4*M4*M  # 2 multiplications
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-2606.43, 2225.51]  [-1907.27, 2779.07]
 [-2700.44, 2863.45]  [-2624.69, 2006.86]

julia> D = M4*M*M4  # 2 multiplications
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-3126.19, 2073.73]  [-2410.01, 2737.56]
 [-2740.19, 2372.8]   [-2575.55, 2563.49]

julia> E = M*M4*M4  # 2 multiplications
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-3077.05, 2643.18]  [-2871.06, 2697.84]
 [-2781.74, 1880.41]  [-2624.69, 2006.86]

julia> F = square(M2*M) * M2*M  # 3 multiplications
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-3310.01, 2898.24]  [-2427.69, 3330.15]
 [-2954.44, 2930.17]  [-2664.27, 2640.69]

julia> G = M2*M * square(M2*M)  # 3 multiplications
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-3432.39, 2898.24]  [-2646.09, 2878.82]
 [-3104.07, 2437.67]  [-2388.99, 2603.88]

julia> H = M2 * square(M2*M) * M  # 3 multiplications
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-3211.14, 2877.44]  [-2416.47, 3160.68]
 [-3121.46, 3105.47]  [-2767.25, 2662.37]

Heuristics

Presumably there is no optimal way to decompose the power, even when ignoring the order of the multiplications. A heuristics could be to decompose into the largest square numbers.

Sharing results

There is another issue: In the above algorithm/heuristics we may compute certain matrix powers several times. It might be better to have some symbolic representation of the decomposition (as we had in the example) and then compute each term only once. However, that sounds more complicated.