Add matrix pencil functions

[mforets]

The following matrix types can be express used to build interval matrix that are only affected by some parameters (ie. parametric interval matrices):

• (A0, A1, p) --> A0 + p A1 --
• (A0, A1, ..., Ak. p) --> A0 + p1 A1 + \lambda_2 A2 + \ldots + pk Ak
• (A0, p, c, A) ---> A0 + p1^{c1} A1 + p2^{c2} A2 + \ldots + pk^{ck} Ak

See #187 for an implementation of the first two.

There is also the closely related concept of matrix pencil, but in that case, the parameter is understood as an indeterminate. I see two options, either:

(a) Use the implementation in #187. Then define a matrix pencil struct with indeterminates. We may or may not compose those in the future.

(b) Define a struct for the (interval) matrix pencil with an indeterminate lambda, which is then instantiated with some p obtaining some other struct which does hold the value of the parameter.

I would go for option (a), mainly because it is an incremental change.

[lucaferranti]

The internal representation in #187 can be ok (whether vector of matrices or matrix of intervalsvectors is better, only time and benchmarks will say), I have a few comments though

1. I think it would be better to separate the domain of the parameters from the parametric structure. The reason for this is that (as you mentioned in zulip) the purpose of this struct would be to model the case where the matrix has some parametric dependencies. I think this model should be independent of the specific value(s) of the interval ranges. If feels a little weird that if I wanted to use the same parametric structure with a different interval range I would have to create a new object. To me, it would make more sense to make those objects callable so that I could do e.g.

A(1..2)
A(0.3)
A([-1..4, 2..3])
A([1.2, 3.4])

etc. which is closer to the mathematical idea of parametric matrix and the notation used in papers. I am ok to accept a "default" cached domain so that A() would evaluate over that, although I find it confusing.

2. Even with the cached default domain, I am not sure I agree with the current implementation of getindex, that is

getindex(M::IntervalMatrixPencil, i::Int) = getindex(M.A0, i) + M.λ * getindex(M.A1, i)

to clarify my previous comment in the PR, I wasn't opposing defining getindex and the array interface, that's definitely a good idea. I was opposing what getindex is currently doing. As commented on that PR, the idea of this struct is to capture the structure A+$\lambda$*B without having to compute the sum and hence losing dependency on parameter(s), while the current definition of getindex is indeed losing that structure and distributing the sum. I think a better solution could be e.g. to return a function, define an AffineExpression struct (as I was doing in the PR in ILA, but doing it properly could require some more work) or resurrect AffineArithmetic.jl (note to myself, update the workflow there) as if we want to use algorithms for parametric linear systems with affine dependency the use of affine arithmetic will become inevitable sooner of later (probably sooner). This would be particularly critical for the case of higher degree polynomials, because in that case you would have dependency problem also inside single elements in the matrix, so you would already overestimate when picking one element. It is true that the general principle of interval arithmetic is "give valid intervals and cope with overestimation", but the whole point of using parametric structures is to fight that. (The third example in the original comment was carefully picked to avoid this problem, but it would also have a very limited application).

3. I am not sure if it was intentional or just delegated to a future PR. But I noticed e.g. that addition between objects of type AffineIntervalMatrix (or whatever the agreed name atm is) is not defined. I don't think we could rely on generic code (which was one selling speech of the current implementation) such as e.g.

+(A::AbstractIntervalMatrix, B::AbstractIntervalMatrix) = map((a, b) -> a + b, A, B)

because given the current implementation of getindex one would lose the parametric structure when accessing the elements in the matrix, hence just computing the sum of two normal interval matrices. One would have to define sum of ParametricIntervalMatrices anyway properly. Something like

+(A::AffineParametricMatrix1, B::AffineParametricMatrix2) = AffineParametricMatrix1(A.A+B.A, A.B + B.B)

(hopefully you get the idea of the previous example, I went a little wild with notation)

4. We should forbid multiplication between two AffineParametricMatrix, because their product would not be a matrix of affine expressions anymore (at least until we have a data structure for higher degree polynomials).

[mforets]

Thanks for your insight.

I think it would be better to separate the domain of the parameters from the parametric structure.

Yes; we can do that through struct composition in subsequent patches. I mean to say, we can add the struct with the indeterminate later, and compose it with the structs in #187.

I think a better solution could be e.g. to return a function,

The current implementation is correct, since the interval parameter is known upon construction. This is connected to the previous comment (to implement a method with a function-like behavior at a later point).

I am not sure if it was intentional or just delegated to a future PR.

Delegated to future PRs. The main motivation of #187 is to start modeling problems with the required structure.