Affine Iterations and Wrapping Effect: Various Approaches

This problem was brought to us through this example of an IIR (Infinite Impulse Response) linear filter in a state-space form:

for $x_{0}=0$ and $x_{1}\in[1,1.1]$. We assume that $u_{n}\in{\mathbf{u}}=[9.95\,,\,10.05]$ for every $n$.

This iteration can also be written as a linear recurrence in $\mathbb{R}^{2}$:

This toy example will be used to illustrate the various approaches mentioned in this paper. The first iterates, obtained using floating-point arithmetic, with random values for $x_{1}\in[1,1.1]$ and each $u_{n}\in[9.95\,,\,10.05]$, are given on the left two columns below.

The system stabilizes around $14$, with variations due to the random values taken by the $u_{n}$. However, the following snippet of Octave code computes the successive iterates using interval arithmetic, using the interval $[1\,,\,1.1]$ for ${\mathbf{x}}_{1}$ and ${\mathbf{u}}=[9.95\,,\,10.05]$ for the $u_{n}$, that is, we replace $u_{n-2}+u_{n-1}+u_{n}$ by $3*{\mathbf{u}}$.

A=[[0 1];[-0.9 1.8]];
xn=[infsup(0,0);infsup(1,1.1)];
b=4.7e-2 * 3.0*[infsup(0,0);infsup(9.95,10.05)];
n=500; for i=1:n, i , xn=A*xn+b, wid(xn(1)), end;

On the right two columns above are the widths of the successive iterates ${\mathbf{x}}_{n}$: the widths of the iterates diverge rapidly to infinity.

The explanation of this phenomenon is the following: the spectral radius of $A$ is strictly less than 1: $\rho(A)\simeq 0.9487<1$, and thus the exact (and, for that matter, floating-point) iterations converge. However, the recurrence satisfied by the widths of the iterates is $\mathrm{wid}({\mathbf{x}}_{n})=1.8*\mathrm{wid}({\mathbf{x}}_{n-1})+0.9*\mathrm{wid}({\mathbf{x}}_{n-2})+4.7.10^{-2}*3*\mathrm{wid}({\mathbf{u}})$, which corresponds to the 2-dimensional iteration $w_{n}=|A|.w_{n-1}+w_{b}$, with $w_{n}=\mathrm{wid}({\mathbf{x}}_{n})$, $|A|$ the matrix whose coefficients are the absolute values of the coefficients of $A$ and $w_{b}=4.7.10^{-2}*3*\mathrm{wid}({\mathbf{u}})$. As the spectral radius of $|A|$ is larger than 1, indeed $\rho(|A|)\simeq 2.208>1$, the iterations diverge.

This phenomenon is a special case of the so-called wrapping effect. Its ubiquity in interval computations has been put in evidence by Lohner in [6].

The wrapping effect is ubiquitous, as defined and developed in [6]. It can be described as the overestimation due to the enclosure of the sought set in a set of a given. simple structure. In our case, this simple structure corresponds to multidimensional intervals or boxes, that is, parallelepipeds with sides parallel to the axes of the coordinate system. When the computation is iterative, and when each iteration produces such an overestimating set that is used as the starting point of the next iteration, the size of the computed set may grow exponentially in the number of iterations, even when the exact solution set remains bounded and small.

Lohner also put in evidence that the affine iteration we study in this paper, namely $x_{n+1}=Ax_{n}+b$, or more generally $x_{n+1}=A_{n}x_{n}+b_{n}$ with $x_{n+1}$, $x_{n}$ and $b_{n}$ vectors in $\mathbb{R}^{d}$ and $A_{n}\in\mathbb{R}^{d\times d}$ for every $n\in\mathbb{N}$, is archetypal. It occurs in many algorithms, and the examples cited in [6] include

• •

  matrix-vector iterations as the ones studied in this paper;

• •

  discrete dynamical systems: ${\mathbf{x}}_{n+1}=f({\mathbf{x}}_{n})$, ${\mathbf{x}}_{0}$ given and $f$ sufficiently smooth;

• •

  continuous dynamical systems (ODEs): $x^{\prime}(t)=g(t,x(t))$, $x(0)=x_{0}$, which is studied through a numerical one step method (or more) of the kind ${\mathbf{x}}_{n+1}={\mathbf{x}}_{n}+h\Phi({\mathbf{x}}_{n},t_{n})+{\mathbf{z}}_{n+1}$;

• •

  difference equations: $a_{0}{\mathbf{z}}_{n}+a_{1}{\mathbf{z}}_{n+1}+\ldots+a_{m}{\mathbf{z}}_{n+m}=b_{n}$ with ${\mathbf{z}}_{1},\ldots{\mathbf{z}}_{m}$ given;

• •

  linear systems with (banded) triangular matrix;

• •

  automatic differentiation.

In this paper, we concentrate on examples similar to the toy example presented above: for every initial value $x_{0}\in{\mathbb{R}}^{d}$, the sequence of iterates $(x_{n})_{n\in{\mathbb{N}}}$ converges to a finite value $x^{*}\in{\mathbb{R}}^{d}$, since $\rho(A)<1$; however, the computations performed using interval arithmetic diverge because their behaviour is dictated by $\rho(|A|)$ which is larger than $1$. We are interested in the iterates computed using interval arithmetic: it is established that these iterates increase in width, however different approaches can be applied to counteract the exponential growth of the width of the iterates. Several of them, some new as in Sections 2.3.2 and 2.3.3, and some already well eestablished as in Section 2.4, will be tried and compared, in terms of the widths of the results and the computational time.

Let $A$ be a $d\times d$ matrix in $\mathbb{R}^{d\times d}$, ${\mathbf{x}}_{0}\in{\mathrm{I}}\!\mathbb{R}^{d}$ be an interval vector (boldface font is used for interval quantities and ${\mathrm{I}}\!\mathbb{R}$ stands for the set of real intervals), $x_{0}$ a vector in $\mathbb{R}^{d}$ with $x_{0}\in{\mathbf{x}}_{0}$, ${\mathbf{b}}\in{\mathrm{I}}\!\mathbb{R}^{d}$ an interval vector, and $b$ a vector in $\mathbb{R}^{d}$ and $b\in{\mathbf{b}}$. In what follows, $n$ denotes the number of iterations.

It is assumed that $\rho(A)<1$ and $\rho(|A|)>1$.

A first goal is to determine the set of all fixed-points of the iteration

for every $x_{0}\in{\mathbf{x}}_{0}$ and every $b\in{\mathbf{b}}$.

It is known that $x_{n}$ can be written as

thus

However, when the vectors $x_{0}$ and $b$ are replaced in the iterative formula by their interval enclosures ${\mathbf{x}}_{0}$ and ${\mathbf{b}}$, one obtains the new interval vector ${\mathbf{x}}_{n+1}$, which is computed as:

Another goal is to determine a tight enclosure for each iterate of this diverging set of intervals.

As mentioned above, the increase in widths of the iterates can be attributed to the use of parallelepipeds with sides parallel to the axes of the coordinate system, and not to the geometry of the transformation. To cure this problem, changes of coordinates will be applied, using an invertible matrix $B$, with $x=By\Leftrightarrow y=B^{-1}x$ and its interval counterpart ${\mathbf{x}}=B{\mathbf{y}}$. This yields the iteration

and its interval counterpart

In what follows, to establish bounds and their proofs, we assume $A$ diagonalizable (this will not necessarily be the case for the experiments) and $A$ can be diagonalized as $A=P^{-1}\Lambda P$ where $\Lambda$ is a diagonal matrix with the eigenvalues $\lambda_{1},\ldots\lambda_{d}$ of $A$ on the diagonal and the columns of $P^{-1}$ are the corresponding eigenvectors.

The iteration considered in this paper corresponds to $x_{n}=A^{n}x_{0}+\sum_{i=0}^{n-1}A^{i}b$. The numerical unstability of computing the matrix power $A^{n}$ and applying it to a vector, is well known: $A^{n}x_{0}$ tends to be aligned with the eigenvector of $A$ associated with the largest (in module) eigenvalue, and the information corresponding to the contribution of the other eigenvectors is lost. To avoid this well-known problem of the power method, we will consider orthogonal changes of coordinates. The choice of the orthogonal matrices is related to $A$, the matrix of the iteration.

We will first consider the QR-factorization of $A$: $A=QR$ with $Q\in{\mathbb{R}}^{d\times d}$ orthogonal, that is, $QQ^{\prime}=Q^{\prime}Q=I$ is the identity matrix and $R\in{\mathbb{R}}^{d\times d}$ is upper triangular.

The other factorization used in this paper is the SVD-factorization of $A$: $A=U\Sigma V^{\prime}$ with $U$, $V$ and $\Sigma$ $\in{\mathbb{R}}^{d\times d}$, where $U$ and $V$ are orthogonal and $\Sigma$ is a diagonal matrix with the singular values $\sigma_{1},\ldots\sigma_{d}$ of $A$ on the diagonal. We also assume that $\sigma_{1}\geq\sigma_{2}\geq\ldots\geq\sigma_{n}\geq 0$. This idea has been sketched but not completely developed by Beaumont in [2].

Mayer and his co-authors have extensively studied the existence of a fixed-point for the iteration studied in this paper. In [7], Mayer and Warnke have thoroughly established formulas for the fixed-point in the case of $\rho(|A|)<1$: this fixed-point is independent of the starting interval ${\mathbf{x}}_{0}$. In [1], Arndt and Mayer have established necessary and sufficient condition on $A$ for a fixed-point to exist, when $\rho(|A|)=1$. In this case, the fixed-point is an interval of nonzero width, that is, a non-degenerate interval. It is well-known that the widths of the iterates diverge when $\rho(|A|)>1$, and thus that no fixed-point exists in this case. Our goal is to study the speed of divergence of the iterates in this case.

The main idea is to use an orthogonal change of coordinates which is related to the matrix of the iteration. As the matrix $A$ is kept constant for all iterations (and this is not the case in the more general approach of Lohner, see Section 2.4), the change of coordinates is also kept constant and given by an orthogonal matrix $B$. The two orthogonal matrices considered in what follows are either $B=Q$ from the QR-factorization of $A$, or $B=U$, resp. $B=V^{\prime}$, from the SVD-factorization of $A$.

A well-known approach is given in Lohner, e.g. in [6] and studied in details by Nedialkov and Jackson in [8]. It is usually presented for the iteration $x_{n+1}=A_{n}x_{n}+b_{n}$, that is when the matrix and the affine term vary at each iteration.

Lohner’s QR method consists in performing the following iteration:

and its interval counterpart is

In the case of a constant – throughout the iterations – matrix $A$, one can recognize Francis’ and Kublanovskaya’s QR-algorithm. Using the convergence of $(R_{n})$ towards the matrix of eigenvalues of $A$ (or towards its Schur form), in [8], Nedialkov and Jackson established the following bounds:

where we recall $A$ diagonalizable: $A=P^{-1}\Lambda P$.

Two aspects are compared: the complexity and the accuracy, that is, the bounds on the widths of the iterates, of each method.

Let us first examine the computational complexity. Let us recall that the QR-factorization, resp. SVD-factorization, of a $d\times d$ matrix has a computational complexity of ${\cal{O}}(d^{3})$. In the algorithms of Sections 2.3.2 and 2.3.3, the factorization of a matrix is performed only once, and not at every iteration: for $n$ iterations, these algorithms thus have complexity ${\cal{O}}(d^{3}+nd^{2})$. In comparison, Lohner’s QR method has complexity ${\cal{O}}(nd^{3})$, which is significantly larger when $d$ is large. In comparison, the cost of the factorization is negligible when the number $n$ of iterations is large.

Let us now compare the accuracy of these different methods, from a theoretical point of view. The bounds we get on the width of the iterate ${\mathbf{x}}_{n}$ are larger than the bounds obtained by Nedialkov and Jackson, as the condition number of the matrix $P$ appears to the $n$-th power in the formula for the QR- and SVD-algorithms, whereas it appears without this $n$-th power in the bound for Lohner’s QR-algorithm. As a condition number is always larger or equal to $1$, this means that the bound for Lohner’s QR-algorithm is tighter than the bounds for the QR- and SVD-algorithms.

After the results on the widths of the iterates of the toy example given in Section 1.1, Section 3.2 presents the computation of each corner of the initial box, to illustrate that it is possible to get tight enclosures, on such a small example.

All algorithms presented in this paper, namely the naive (or brute-force) application of the iteration, the QR-algorithm of Section 2.3.2, the two versions of the SVD-algorithm of Section 2.3.3, and Lohner’s QR algorithm given in Section 2.4 have been implemented in Octave using Heimlich’ interval package [3], then in Matlab using Rump’s Intlab package [10]. Two other methods have been implemented and compared. The first technique [9] consists in the determination of $k$ such that $\rho(|A^{k}|)<1$, then it computes only one iterate every $k$ step, in other words it computes

this iteration converges even when interval arithmetic is employed. The other technique is the use of affine arithmetic, as advocated by Rump in a private communication. In the experimental results presented below, each technique is associated to a color:

The factorizations use only the basic QR and SVD factorizations available in Matlab, but neither the pivoted QR recommended by Lohner in [6] nor more elaborate versions presented by Higham in [4].

Sections 3.3 and 3.4 contain the evolution of the radii of the iterates computed by these different techniques, for two matrix dimensions: $10\times 10$ and $100\times 100$. The $y$-axis for the radii uses a logarithmic scale. For both dimensions, four kinds of matrices $A$ have been used for the experiments. On the one hand, matrices which are well-conditioned (with a condition number of order $10^{2}$) and ill-conditioned (with a condition number of order $10^{10}$) have been generated. On the other hand, the scaling of the matrices varies: matrices which are well-scaled and matrices which are ill-scaled, with the order of magnitude of their coefficients varying between $1$ and $10^{10}$. These are only orders of magnitudes, as the matrices, originally generated by a call to Matlab’s randsvd were then added to a multiple of the identity matrix and multiplied by a constant, in order to satisfy both $\rho(A)<1$ and $\rho(|A|)>1$. It can also be noted that degrading the scaling of the matrix also degrades its condition number; in other words, a “well-conditioned ill-scaled” matrix has a much worse condition number than a “well-conditioned well-scaled” matrix, even if the required condition numbers, in the call to randsvd, are initially the same.

All experiments have been performed on a 2.7 GHz Quad-Core Intel Core i7 with 16GB RAM. Timings are averaged over 100 executions, except for affine arithmetic where at most 10 executions were performed.

First, the toy example presented in Section 1.1 is considered. As the iteration is affine, one can compute separately the images of the endpoints of the initial vector, to get the endpoints of the successive iterates. That is, we compute separately

for $x_{0}=0$ and $x_{1}=1$ and for $x_{0}=0$ and $x_{1}=1.1$. However, we use the interval vector ${\mathbf{u}}=[9.95,10.05]$ in the iteration. The convex hull of the 10 first iterates are represented on the left part of Figure 1. It is obvious that the width of the successive iterates grow rapidly.

Then we compute separately

for $x_{0}=0$ and $x_{1}=1$ and for $x_{0}=0$ and $x_{1}=1.1$, and for $u=9.95$ and $u=10.05$. The convex hull of the 10 first iterates are represented on the right part of Figure 1. In this case, the width of the successive iterates remain small, of the order of magnitude of $1\%$ of the midpoint of the interval.

In this toy example, the iterations had to be performed 4 times, that is, once for each corner of the initial values, in order to get a tight enclosure. This can clearly not be generalized to high dimensions, as the number of corners grows as $4^{d}$ with the dimension $d$ of the problem.

Figure 2 gives the radii (in logarithmic scale) for the successive iterates computed by the methods detailed above. When the number of iterations is large (visually, above 30 or 40 iterations), the iterates computed by all methods presented in Section 2 diverge rapidly, as can be seen on the plots on the left. When one concentrates on the first iterations, the behaviours compare differently. One can also note that unscaling the matrix $A$ speeds the divergence, for all methods. On the contrary, the $k$-step method and the use of affine arithmetic preserve the convergence of the iterates.

The timings in seconds are given below:

The methods presented in Sections 2.3.2, 2.3.3 and 2.4 all exhibit similar execution times. The naive method performs less operations and is thus faster. The $k$-th step method is fast as well, the variations in its execution time are due to the preprocessing, that is to the determination of the power $k$ such that $\rho(|A^{k}|)<1$: the execution time is larger when $k$ is larger. With this method, the convergence is good. The use of affine arithmetic significantly slows down the computations, however the iterates converge.

Figure 3 gives the radii (in logarithmic scale) for the successive iterates computed by the methods detailed in Section 3.1. When the number of iterations is large (visually, above 40 or 50 iterations), the iterates computed by all methods, except the $k$-step method, diverge rapidly, as can be seen on the plots on the left. Again, when one concentrates on the first iterations, the behaviours compare differently.

The timings in seconds are given below:

The comments on the timings apply again, with the exception of the use of affine arithmetic, which is still much slower but does not manage any more to preserve the convergence very long.

One can note that the $k$-step method, that is the method that resorts to a convergent interval iteration, performs very well at a moderate computation cost. Even the preprocessing time to determine the value of $k$ has a negligible cost.

This method is a totally ad hoc approach for this problem and cannot be generalized. However, in the framework of filters and control theory, it has a physical meaning: the divergence of the iterations can be attributed to a sampling time which is too small to allow variations to be observed. Multiplying the sampling time by $k$ means sampling less frequently (by a factor $k$) and thus being able to measure the evolution of the observed quantities.

The use of affine arithmetic, on the contrary, is a very general method and it exhibits a very good accuracy, even if it eventually diverges (see the experiments with the $100\times 100$ matrices in Section 3.4). The counterpart is the execution time, which is at least a thousand times larger than for the other methods. This is not an issue for the experiments presented here, as the time is of order of magnitude of a minute.

The methods based on the QR or SVD factorizations of the matrix $A$ were developed with geometric principles in mind. For the QR-algorithm, the idea was to align the current box with the directions that are preserved by the product by $A$, with a tradeoff between aligning the box along the eigenvectors and preserving an orthonormal system of coordinates, hence the choice of $Q$. For the SVD-algorithm, the idea was to align the box along the direction which gets the maximal elongation, that is along the singular vector corresponding to the largest singular value.

In both cases, the benefit of these geometric transformations is mitigated with the overestimation implied by extra computations, and there is either no clear benefit for the QR-based approach, or a delicate balance for the SVD-based approach. The SVD-algorithm is interesting when the matrix is ill-scaled, and particularly for the first iterations.

The methods of choice remain either the naive approach, when the matrix $A$ is well-conditioned and well-scaled, or Lohner’s QR method when the matrix is ill-conditioned. Surprisingly, the overhead of Lohner’s QR method, in terms of computational time, is not as large as the formula for its complexity implies.

Our general recommendation is thus:

• •

  to preprocess the matrix $A$ in order to scale it;

• •

  then to execute in parallel the naive approach and Lohner’s QR approach, in order to converge reasonably well for any condition number of $A$.

Affine arithmetic is a solution of choice when other solutions fail and when the analysis and developing time is a scarce resource.