Numerical computation of formal solutions to interval linear systems of equations

TL;DR

This paper develops numerical methods using Kaucher interval arithmetic to compute formal solutions of interval linear systems, addressing convergence and implementation challenges.

Contribution

It introduces two stationary iterative approaches based on matrix splitting for formal solutions in Kaucher arithmetic, comparing their effectiveness.

Findings

Two iterative methods successfully compute formal solutions.

Convergence conditions are established for the proposed methods.

Comparison shows advantages over existing approaches.

Abstract

The work is devoted to the development of numerical methods for computing "formal solutions" of interval systems of linear algebraic equations. These solutions are found in Kaucher interval arithmetic, which extends and completes the classical interval arithmetic algebraically. The need to solve these problems naturally arises, for example, in inner and outer estimation of various solution sets to interval linear systems of equations. The work develops two approaches to the construction of stationary iterative methods for computing the formal solutions that are based on splitting the matrix of the system. We consider their convergence and implementation issues, compare with the other approaches to computing formal solutions.

Tables

Table 1. Table 1: Multiplication in Kaucher complete interval arithmetic

Table 2. Table 2: Algorithm TrnSplit for computing formal solutions,

Input

An interval linear algebraic $n\times n$-system $\text{A}x=\text{b}$.

A triangular splitting of the matrix A of the system

I to interval matrices $\text{G}=\left({\text{g}}_{ij}\right)$ and $\text{H}=\left({\text{h}}_{ij}\right)$.

A specified accuracy $ϵ$.

Output

An approximation to formal solution of the system $\text{A}x=\text{b}$.

Algorithm

$q\leftarrow +\infty$;       assign a starting value to the vector x;       DO WHILE   ( $q\ge ϵ$ )       ${\text{p}}_1\leftarrow {\text{b}}_1$;       DO FOR  $i=2$  TO  $n$       ${\text{p}}_i\leftarrow {\text{b}}_i\ominus \sum _{j=1}^{i-1}{\text{h}}_{ij}{\text{x}}_j$       END DO       ${\tilde{\text{x}}}_n\leftarrow {\text{p}}_n\oslash {\text{g}}_{nn}$;       DO FOR  $i=n-1$  TO  $1$   STEP  $\left(-1\right)$       ${\tilde{\text{x}}}_i\leftarrow \left({\text{p}}_i\ominus \sum _{j=i+1}^n{\text{g}}_{ij}{\tilde{\text{x}}}_j\right)\oslash {\text{g}}_{ii}$       END DO       $q\leftarrow$ distance between the vectors x and $\tilde{\text{x}}$;       $\text{x}\leftarrow \tilde{\text{x}}$;       END DO

Equations

$$\left\{\begin{array}[]{ccccccccc}\textbf{{a}}_{11}x_{1}&+&\textbf{{a}}_{12}x_{2}&+&\ldots&+&\textbf{{a}}_{1n}x_{n}&=&\textbf{{b}}_{1},\\[5.69054pt] \textbf{{a}}_{21}x_{1}&+&\textbf{{a}}_{22}x_{2}&+&\ldots&+&\textbf{{a}}_{2n}x_{n}&=&\textbf{{b}}_{2},\\[5.69054pt] &\vdots&&&\ddots&&&\vdots&\\[5.69054pt] \textbf{{a}}_{n1}x_{1}&+&\textbf{{a}}_{n2}x_{2}&+&\ldots&+&\textbf{{a}}_{nn}x_{n}&=&\textbf{{b}}_{n},\end{array}\right.$$

$$\text{A}x=\text{b},$$

$$\text{a}\star \text{b}=\{a\star b\mid a\in \text{a},b\in \text{b}\}.$$

$$\text{a}+\text{b}=\left[\underline{\text{a}}+\underline{\text{b}},\overline{\text{a}}+\overline{\text{b}}\right],$$

$$\text{a}-\text{b}=\left[\underline{\text{a}}-\overline{\text{b}},\overline{\text{a}}-\underline{\text{b}}\right],$$

$$\text{a}\cdot \text{b}=[\min \{\underline{\text{a}}\underline{\text{b}},\underline{\text{a}}\overline{\text{b}},\overline{\text{a}}\underline{\text{b}},\overline{\text{a}}\overline{\text{b}}\},\max \{\underline{\text{a}}\underline{\text{b}},\underline{\text{a}}\overline{\text{b}},\overline{\text{a}}\underline{\text{b}},\overline{\text{a}}\overline{\text{b}}\}],$$

$$\text{a}/\text{b}=\text{a}\cdot \left[1/\overline{\text{b}},1/\underline{\text{b}}\right]\text{mbox}for\text{b}\not\ni 0.$$

$$\text{a}+x=\text{b},\text{a}\cdot x=\text{b}$$

$$\text{a}\wedge \text{b}:=$$

$$\subseteq$$

$$\text{a}\vee \text{b}:=$$

$$\subseteq$$

$$\mathrm{dual}\text{a}:=[\overline{\text{a}},\underline{\text{a}}].$$

$${\mathrm{pro}\;}\textbf{{a}}:=\begin{cases}\ \textbf{{a}},&\mbox{ if $\textbf{{a}}$ is proper, }\\[2.84526pt] \ {\mathrm{dual}\;}\textbf{{a}},&\mbox{ otherwise. }\end{cases}$$

$$\text{a}\subseteq \text{b}⟺\underline{\text{a}}\ge \underline{\text{b}}\text{mbox}and\overline{\text{a}}\le \overline{\text{b}}.$$

$${\mathrm{sgn}\;}\textbf{{a}}=\left\{\ \begin{array}[]{cl}+\ ,&\mbox{ if }\;\textbf{{a}}\geq 0,\\[2.84526pt] -\ ,&\mbox{ if }\;\textbf{{a}}\leq 0,\\[2.84526pt] \mbox{undefined},&\mbox{ if $\,0\,$ is in the interior of ${\mathrm{pro}\;}\textbf{{a}}$}.\end{array}\right.$$

$$\text{a}+\text{b}:=\left[\underline{\text{a}}+\underline{\text{b}},\overline{\text{a}}+\overline{\text{b}}\right].$$

$$\mathrm{opp}\text{a}:=[-\underline{\text{a}},-\overline{\text{a}}].$$

$$\text{a}\ominus \text{b}:=\text{a}+\mathrm{opp}\text{b}=\left[\underline{\text{a}}-\underline{\text{b}},\overline{\text{a}}-\overline{\text{b}}\right].$$

$$\mathrm{dual}(\text{a}+\text{b})=\mathrm{dual}\text{a}+\mathrm{dual}\text{b}.$$

$$\mu\cdot\textbf{{a}}\ :=\ \left\{\ \begin{array}[]{ll}[\,\mu\,\underline{\textbf{{a}}},\mu\,\overline{\textbf{{a}}}\,],&\mbox{ if }\;\mu\geq 0,\\[2.84526pt] [\,\mu\,\overline{\textbf{{a}}},\mu\,\underline{\textbf{{a}}}\,],&\mbox{ otherwise.}\end{array}\right.$$

$$\mathrm{dual}(\text{a}\cdot \text{b})=\mathrm{dual}\text{a}\cdot \mathrm{dual}\text{b}.$$

$$\displaystyle\textbf{{a}}\star\textbf{{b}}\$$

$$\displaystyle=\ \raisebox{-0.86108pt}{$\Bigl{[}$}\;\min_{a\in\textbf{{a}}}\,\min_{b\in\textbf{{b}}}\,(a\star b),\ \max_{a\in\textbf{{a}}}\,\max_{b\in\textbf{{b}}}\,(a\star b)\;\raisebox{-0.86108pt}{$\Bigr{]}$}=\ \bigvee_{a\in\textbf{{a}}}\ \bigvee_{b\in\textbf{{b}}}\ (a\star b),$$

$$\text{a}\star \text{b}=\underset{a\in \mathrm{pro}\text{a}}{\bigwedge }\underset{b\in \mathrm{pro}\text{b}}{\bigwedge }(a\star b),$$

$$\text{a}\star \text{b}=\underset{a\in \text{a}}{\bigvee }\underset{b\in \mathrm{pro}\text{b}}{\bigwedge }(a\star b),$$

$$\text{a}\star \text{b}=\underset{a\in \mathrm{pro}\text{a}}{\bigwedge }\underset{b\in \text{b}}{\bigvee }(a\star b),$$

$$\begin{array}[t]{c}\\[-15.0pt] \mbox{\begin{picture}(5.0,5.0)(-1.0,0.0) \put(0.0,5.0){\line(0,-1){5.0}} \put(3.0,5.0){\line(0,-1){5.0}} \put(0.2,0.0){\line(3,5){3.0}}\put(0.1,0.0){\line(3,5){3.0}} \put(0.0,0.0){\line(3,5){3.0}} \put(0.1,0.0){\line(3,5){3.0}} \put(0.2,0.0){\line(3,5){3.0}}\end{picture}}^{\scriptstyle\textbf{{a}}}\\[-7.0pt] {\scriptstyle x}\end{array}:=\left\{\ \begin{array}[]{cl}\ \displaystyle\bigvee_{x\in\textbf{{a}}},&\mbox{ if $\textbf{{a}}$ is proper,}\\[17.07164pt] \ \displaystyle\bigwedge_{x\in{\mathrm{dual}\;}\textbf{{a}}},&\mbox{ if $\textbf{{a}}$ is improper.}\end{array}\right.$$

$$\textbf{{a}}=\begin{array}[t]{c}\\[-15.0pt] \mbox{\begin{picture}(5.0,5.0)(-1.0,0.0) \put(0.0,5.0){\line(0,-1){5.0}} \put(3.0,5.0){\line(0,-1){5.0}} \put(0.2,0.0){\line(3,5){3.0}}\put(0.1,0.0){\line(3,5){3.0}} \put(0.0,0.0){\line(3,5){3.0}} \put(0.1,0.0){\line(3,5){3.0}} \put(0.2,0.0){\line(3,5){3.0}}\end{picture}}^{\scriptstyle\textbf{{a}}}\\[-7.0pt] {\scriptstyle x}\end{array}x$$

$$\textbf{{a}}\star\textbf{{b}}\ =\ \begin{array}[t]{c}\\[-15.0pt] \mbox{\begin{picture}(5.0,5.0)(-1.0,0.0) \put(0.0,5.0){\line(0,-1){5.0}} \put(3.0,5.0){\line(0,-1){5.0}} \put(0.2,0.0){\line(3,5){3.0}}\put(0.1,0.0){\line(3,5){3.0}} \put(0.0,0.0){\line(3,5){3.0}} \put(0.1,0.0){\line(3,5){3.0}} \put(0.2,0.0){\line(3,5){3.0}}\end{picture}}^{\scriptstyle\textbf{{a}}}\\[-7.0pt] {\scriptstyle a}\end{array}\begin{array}[t]{c}\\[-15.0pt] \mbox{\begin{picture}(5.0,5.0)(-1.0,0.0) \put(0.0,5.0){\line(0,-1){5.0}} \put(3.0,5.0){\line(0,-1){5.0}} \put(0.2,0.0){\line(3,5){3.0}}\put(0.1,0.0){\line(3,5){3.0}} \put(0.0,0.0){\line(3,5){3.0}} \put(0.1,0.0){\line(3,5){3.0}} \put(0.2,0.0){\line(3,5){3.0}}\end{picture}}^{\scriptstyle\textbf{{b}}}\\[-7.0pt] {\scriptstyle b}\end{array}\ (a\star b).$$

$$\mathcal{P}$$

$$\mathcal{Z}$$

$$-\mathcal{P}$$

$$\mathrm{dual}\mathcal{Z}$$

$$\text{ab}=\text{ac}\Rightarrow \text{b}=\text{c}.$$

Full text

1 Introduction

The main object we study in this article is interval systems of linear algebraic equations having the form

[TABLE]

with intervals $\textbf{{a}}_{ij},\textbf{{b}}_{i}$, which will be considered as elements of Kaucher complete interval arithmetic $\mathbb{KR}$ [16]. It is convenient to denote such systems as

[TABLE]

where $\textbf{{A}}=(\,\textbf{{a}}_{ij})$ is an interval $n\times n$-matrix, $\textbf{{b}}=(\,\textbf{{b}}_{i})$ is an interval $n$-vector.

For interval equations and systems of equations, various definitions of solutions and solution sets exist. Additionally, practice requires that these solution sets be evaluated in a wide variety of ways. As a consequence, there are many problem statements associated with interval equations and systems of equations. In our paper, we are going to concentrate on computing the so-called formal solutions to interval linear systems of the form (1)–(2).

Definition 1* An interval vector is called formal solution to an interval equation or a system of equations, if substituting it into the equation or equation system and execution of all operations in interval arithmetic result in a valid equality.*

Due to a number of reasons, it is useful to find formal solutions to interval systems of equations in Kaucher interval arithmetic $\mathbb{KR}$ rather than in classical interval arithmetic $\mathbb{IR}$. $\mathbb{KR}$ is an algebraic and order completion of $\mathbb{IR}$.

The problem we consider is not new. It was first noted in the paper by S. Berti [2] with respect to one interval quadratic equation. Then H. Ratschek and W. Sauer [32] studied such solutions for a single interval linear equation, and they used the term algebraic solution. K. Nickel [29] considered formal solutions to interval linear systems of equations in complex interval arithmetics. However, formal (algebraic) solutions of interval equations and systems of equations have been studied only in about a dozen papers for the last decades when interval analysis rapidly developed.

The need to consider formal solutions arises, in particular, in the interval analogue of the method of undetermined coefficients, as well as in modeling linear static systems under interval uncertainty of their parameters, as a consequence of the so-called formal (algebraic) approach to their analysis [36, 37, 39]. We remind that this approach replaces the original problem of estimating the internal states of the system to the problem of finding algebraic solutions of an auxiliary equation in Kaucher arithmetic $\mathbb{KR}$. In this work, we do not focus on the derivation of Kaucher interval arithmetic and discussion of its remarkable properties. The interested reader should refer to the original works of [9, 10, 16], or to the summary in [20, 36, 37, 38, 39].

Everywhere below, we will denote the intervals and interval objects (vectors, matrices) in bold (for example, A, B, C, …, x, y, z), $\underline{\textbf{{x}}}$, $\overline{\textbf{{x}}}$ are understood as the lower (left) and upper (right) endpoints of an interval x, other symbols also follow the informal international standard [18]. Some results of this work were previously published in [36].

2 Kaucher interval arithmetic

The classical interval arithmetic $\mathbb{IR}$ is known to be an algebraic system formed by intervals of the real axis $\mathbb{R}$ with the operations between them defined “by representatives”, i. e. according to the following fundamental principle

[TABLE]

The constructive reformulation of the above relation for separate arithmetic operations looks as follows:

[TABLE]

The algebraic properties of the classical interval arithmetic $\mathbb{IR}$ are much more poor than those of the usual number systems, the ring $\mathbb{Z}$ of integers, the fields of rational numbers $\mathbb{Q}$ and real numbers $\mathbb{R}$, because

• •

all intervals with nonzero width, i. e. most of the elements of $\mathbb{IR}$ do not have inverse elements with respect to the operations (4)–(7),

• •

the arithmetic operations (4)–(7) are related to each other by a weak sub-distributivity relation (see [27, 28, 40]), and the full distributivity of multiplication (and division) with respect to addition and subtraction does not take place.

As a consequence, first, in $\mathbb{IR}$ elementary equations with respect to the unknown variable $x$

[TABLE]

and their like do not always have solutions. Secondly, the technique of symbolic transformations in the classical interval arithmetic $\mathbb{IR}$ is quite poor. We cannot even transfer terms from one part of the equation to another, and the lack of distributivity makes it impossible to reduce such terms.

In addition, the ordinal properties of the classical interval arithmetic with respect to inclusion ordering “$\subseteq$” are unsatisfactory. In partially ordered sets, the possibility of taking, for any two elements, their lower bound “$\wedge$” and the upper bound “$\vee$” with respect to the order in question plays a huge role. In $\mathbb{IR}$, the corresponding operations are

[TABLE]

Since the first of these operations is not always fulfilled, then $\mathbb{IR}$ is, in a sense, “not closed” with respect to the inclusion order.111If $\textbf{{a}},\textbf{{b}}$ are usual one-dimensional intervals with non-empty intersection, then $\textbf{{a}}\wedge\textbf{{b}}$ and $\textbf{{a}}\vee\textbf{{b}}$ coincide with $\textbf{{a}}\cap\textbf{{b}}$ and $\textbf{{a}}\cup\textbf{{b}}$ respectively. But this is not valid in general. For example, $[1,2]\wedge[3,4]$ is not defined.

The elements of the complete interval arithmetic $\mathbb{KR}$ are pairs of real numbers $[\,\eta,\,\vartheta\,]$, not necessarily connected by the relation $\,\eta\leq\vartheta$. Thus, $\mathbb{KR}$ is obtained by appending improper intervals $[\,\eta,\vartheta\,]$, $\eta>\vartheta$ to the set $\mathbb{IR}=\{\,[\,\eta,\vartheta\,]\mid\eta,\,\vartheta\in\mathbb{R},\;\eta\leq\vartheta\,\}$ of proper intervals and real numbers (identified with degenerate intervals of zero width). The elements of Kaucher arithmetic and the more complex objects formed from them (vectors, matrices) will be highlighted in bold, similar to usual intervals. Moreover, if $\textbf{{a}}=[\,\eta,\,\vartheta\,]$, then $\eta$ is called the left (or lower) end of the interval a and is denoted by $\underline{\textbf{{a}}}$ or $\inf\textbf{{a}}$, and $\vartheta$ is called the right (or upper) endpoint of the interval a and is denoted $\overline{\textbf{{a}}}$ or $\sup{a}$. The interval a is called balanced if $\underline{\textbf{{a}}}\ =\ -\overline{\textbf{{a}}}$.

Without trying to replace the rigorous mathematical constructions performed by E. Kaucher, we will consider the description of the arithmetic $\mathbb{KR}$ along with informal motivations that help to understand the meaning of various constructions.

Proper and improper intervals, the two halves of $\mathbb{KR}$, transform into each other as a result of the dualization mapping $\,{\mathrm{dual}\;}:\mathbb{KR}\to\mathbb{KR}$ swapping (turning over) the endpoints of the interval, i. e. such that

[TABLE]

Proper projection of an interval a is the value

[TABLE]

Similar to the classical interval arithmetic $\mathbb{IR}$, the “inclusion” of one interval to another is defined in $\mathbb{KR}$ as follows:

[TABLE]

For example, $[3,1]\subseteq[2,2]=2\in\mathbb{R}$. As a consequence, the operations of taking the minimum (8) and the maximum (9) keep their definitions unchanged in $\mathbb{KR}$, but now they are always possible due to the presence of improper intervals. In particular, $[1,2]\wedge[3,4]=[3,2]$. Thus, the extension of $\mathbb{IR}$ to $\mathbb{KR}$ makes the set of intervals a lattice, and even a conditionally complete lattice with respect to the inclusion ordering (10).222A conditionally complete lattice is a partially ordered set in which every non-empty bounded subset has exact upper and lower bounds [3]. This is more than just a lattice, but less than a full lattice, where minima and maxima can be taken for any families of elements.

In addition to the set-theoretic inclusion on the set of intervals $\mathbb{KR}$, there is another partial ordering, which naturally generalizes the linear order “$\leq$” on the real axis:

Definition 2* For the intervals a, $\textbf{{b}}\in\mathbb{KR}$, we say that a does not exceed b and write “$\,\textbf{{a}}\leq\textbf{{b}}$” if and only if $\,\underline{\textbf{{a}}}\leq\underline{\textbf{{b}}}\,$ and $\,\overline{\textbf{{a}}}\leq\overline{\textbf{{b}}}$.*

The interval is called nonnegative, i. e. “$\,\geq 0$” if both its endpoints are nonnegative. The interval is nonpositive, i. e. “$\,\leq 0$” if both its endpoints are not positive.

For example, $[1,2]\leq[3,2]$, and both compared intervals $[1,2]$ and $[3,2]$ are non-negative.

It is useful to introduce the concept of the sign of an interval, which we define as

[TABLE]

The zero, i. e. the zero interval $[0,0]$, may have any sign.

The semigroup of all proper intervals with the operation of addition is fairly simple: the addition of intervals is divided into independent operations of addition of the left and right endpoints of the operands. As a consequence, the extension of addition from $\mathbb{IR}$ to $\mathbb{KR}$ is easy, and it is defined in $\mathbb{KR}$ in exactly the same way as in classical interval arithmetic:

[TABLE]

But now it follows from the existence of improper intervals that each element a of $\mathbb{KR}$ has a unique inverse with respect to addition (also called opposite), denoted by “${\mathrm{opp}\;}\textbf{{a}}$”, and the equality $\,\textbf{{a}}+{\mathrm{opp}\;}\textbf{{a}}=0\,$ implies that

[TABLE]

With respect to addition, the arithmetic $\mathbb{KR}$ is thus a commutative group that is isomorphic to the additive group of the standard linear space $\mathbb{R}^{2}$. For brevity, we denote by “$\ominus$” an operation that is the inverse of addition, and it will be called internal subtraction in $\mathbb{KR}$ (or algebraic subtraction). Then

[TABLE]

It is easy to verify that, for the addition in $\mathbb{KR}$, the following relation is valid

[TABLE]

It means, in particular, that the addition of improper intervals in $\mathbb{KR}$ is a “mirror image” of the addition for proper intervals in $\mathbb{IR}$. This is consistent with the status of improper intervals, so it seems necessary to preserve this property when defining other arithmetic operations in $\mathbb{KR}$.

We start extending the multiplication to the complete interval arithmetic $\mathbb{KR}$ with the simplest situation, i. e. multiplication by real numbers. It is natural to define it in exactly the same way as in classical interval arithmetic:

[TABLE]

Further, for nonnegative proper intervals, multiplication in $\mathbb{IR}$ is performed “through the endpoints”, quite similar to addition. Therefore, it makes sense to define, in $\mathbb{KR}$, the multiplication of non-negative intervals, which are not necessarily proper, according to the same “endpoint formulas”, i. e. as $\textbf{{a}}\cdot\textbf{{b}}=[\,\underline{\textbf{{a}}}\,\underline{\textbf{{b}}},\overline{\textbf{{a}}}\,\overline{\textbf{{b}}}\,]$. In fact, in this way, we embed a multiplicative semigroup of positive intervals into a group. If, in the product $\textbf{{a}}\cdot\textbf{{b}}$, one of the factors is a nonpositive interval, then we can make it non-negative through multiplying by $(-1)$ and using formula (13), thus coming to the previous case.

Note that in all these cases the property analogous to (12) holds:

[TABLE]

How to extend the definition of multiplication to the entire set $\mathbb{KR}$? The ability of algebra to do this has already been exhausted, and we need to involve considerations concerning the inclusion ordering in $\mathbb{KR}$ and the related properties of arithmetic operations.

Using the maximum with respect to inclusion (9), the fundamental property (3), which defines the operations of classical interval arithmetic, can be rewritten in the following equivalent form:

[TABLE]

where $\star\in\{\,+,-,\cdot\,,/\,\}$. It is easily seen that the addition extended to the entire set of proper and improper intervals $\mathbb{KR}$, as well as the multiplication defined for intervals that do not contain zero and are not contained in zero can be represented in a similar way through the operations (8) and (9) of taking minimum and maximum with respect to inclusion. If both operands a and b are improper, then

[TABLE]

where $\star\in\{+,\cdot\}$. The lower arguments of the operations “$\wedge$” should have proper projections ${\mathrm{pro}\;}\textbf{{a}}$ and ${\mathrm{pro}\;}\textbf{{b}}$, since improper intervals themselves are contained in points due to the definition of inclusion in $\mathbb{KR}$.

The property (15) is evidently obtained from (14) with the use of dualization. But the next properties of the operation $\star\in\{\,+,\cdot\,\}$ are not obvious:

[TABLE]

if a is proper and b is improper, and

[TABLE]

if a is improper and b is proper. Their validity can be verified by direct verification.

All the above written formulas, starting with (14), can be combined into one as follows. We introduce the so-called conditional operation of taking the extremum with respect to inclusion:

[TABLE]

This is an operation that depends on the interval parameter a standing as its upper index. The operation is either maximum or minimum with respect to the inclusion “$\subseteq$”, depending on whether a is proper or improper. This extremum is taken over all $x$ from the proper projection of the interval a. Note that any interval $\textbf{{a}}\in\mathbb{KR}$ can be represented as

[TABLE]

(instead of $x\,$, any letter can be used in the formula). Anyway, for $\star\in\{+,\cdot\}$, the following relation is valid:

[TABLE]

This representation, first obtained in [14], expresses the relationship between the result of the interval operation $\textbf{{a}}\star\textbf{{b}}\,$ and the results of the point operations $a\star b$ for $a\in{\mathrm{pro}\;}\textbf{{a}}$ and $b\in{\mathrm{pro}\;}\textbf{{b}}$. It can be taken as a basis for the definition of arithmetic operations in the complete interval arithmetic $\mathbb{KR}$.

It is not hard to derive, from (17), the monotonicity of the interval arithmetic operations with respect to inclusion.

In order to write out explicit formulas for multiplication in complete interval arithmetic, we select the following subsets in $\mathbb{KR}$:

[TABLE]

Overall, $\mathbb{KR}=\mathcal{P}\,\cup\,\mathcal{Z}\,\cup\,(-\mathcal{P})\,\cup\,({\mathrm{dual}\;}\mathcal{Z})$. Then the multiplication in Kaucher interval arithmetic can be described by Table 1 [16], the cells of which are obtained as a result of detailed writing out the particular cases of applying formula (17) and our previous results. A remarkable fact is that this table is the supplement of a similar table for multiplication in classical interval arithmetic with one more row and one more column that correspond to the case of operands from the set ${\mathrm{dual}\;}\mathcal{Z}$.

As we see, the multiplication in Kaucher arithmetic admits non-trivial zero divisors. For example, $[\,-1,\,2\,]\cdot[\,5,\,-3\,]=0$. The interval multiplication in Kaucher arithmetic turns out to be commutative and associative [9, 15, 16], but the multiplication group in $\mathbb{KR}$ is formed only by intervals a for which $\underline{\textbf{{a}}}\,\overline{\textbf{{a}}}>0$ (or, otherwise, $0\not\in{\mathrm{pro}\;}\textbf{{a}}$), because no any wider subset of $\mathbb{KR}$ satisfies the so-called “cancellation law”

[TABLE]

This is the algebraic condition that a semigroup can be embedded into a group.

Therefore, for any interval a of $\mathbb{KR}$ that does not contain zero and is not contained in zero itself, there is a single inverse element with respect to multiplication, which we will denote by “${\mathrm{inv}\,}\textbf{{a}}$”. From the equality $\,\textbf{{a}}\cdot{\mathrm{inv}\,}\textbf{{a}}=1\,$, it follows that

[TABLE]

For brevity, we will denote the inverse operation of the multiplication, the so-called internal (algebraic) division in $\mathbb{KR}$, by “$\oslash$”, so that

[TABLE]

The above table of explicit formulas for multiplication in the complete interval arithmetic is convenient for computer implementation, but rather cumbersome and barely foreseeable. In some cases, it makes sense to resort to other formulas for interval multiplication in $\mathbb{KR}$, which were proposed by A.V. Lakeyev in [21, 23]. Recall the following definition (see, for example, [3, 5]):

Definition 3*

For a real number $a$, the values*

[TABLE]

are called positive part and negative part of $a$ respectively.

Then $\,a=a^{+}-a^{-}$ and $|a|=a^{+}+a^{-}$.

Proposition 1* (Lakeyev formulas)  For any intervals  a, $\textbf{{b}}\in\mathbb{KR}$, there holds*

[TABLE]

If one of the intervals is a, b is proper, then

[TABLE]

This formula is not simplified if we additionally know that both intervals a, b are proper.

If one of the intervals a, b is proper and the other is improper, then

[TABLE]

A**

The advantage of the Lakeyev formulas is their global character. They give a single and uniform expression for the interval product $\textbf{{a}}\cdot\textbf{{b}}$ over the entire domain of a and b, whereas the representation via Table. 1 has a piecewise character. This is inconvenient in the study of interval functions “as a whole”, in particular, in the study of differentiability and its analogues, in the calculation and evaluation of generalized derivatives, etc.

Subtraction and division in arithmetic $\mathbb{KR}$ are defined in the same way as in classical interval arithmetic:

[TABLE]

Similar to its classical predecessors, all operations of the complete interval arithmetic are inclusion monotone, i. e. monotone with respect to the partial order (10): if $\textbf{{a}},\textbf{{a}}^{\prime},\textbf{{b}},\textbf{{b}}^{\prime}\in\mathbb{KR}$, then

[TABLE]

for any arithmetic operation $\star\in\{\;+,\,-,\,\cdot\,,\,/\;\}$. This follows from their definition according to formula (17).

The relationship of addition and multiplication in Kaucher arithmetic is expressed by the following inclusions:

[TABLE]

These inclusions turn to exact equalities when, in particular, a squeezes to a point, that is, $\textbf{{a}}=a\in\mathbb{R}$:

[TABLE]

Another important case of distributivity is the case when the signs of the intervals b and c coincide with each other:

[TABLE]

E. Gardeñes et al. introduced in [9], for a complete description of all cases of distributivity, the concept of distributive areas defined by an interval a. Membership of the operands in these distributivity areas leads to equalities instead of inclusions (21)–(22). Later, S. Markov and his co-workers, as a classification of various particular cases of the distributivity of addition with respect to multiplication in $\mathbb{KR}$, proposed a “generalized distribution law” [6, 26, 31], covering a large number of various situations. Of the variety of cases considered in these articles, we will further need the following relation:

[TABLE]

if the intervals  b, c and $\textbf{{b}}+\textbf{{c}}\,$ have definite signs and $\,{\mathrm{sgn}\;}\textbf{{b}}=-{\mathrm{sgn}\;}\textbf{{c}}={\mathrm{sgn}\;}(\textbf{{b}}+\textbf{{c}})$.

3 Theory

In this section, we consider the basic theoretical facts concerning interval linear systems of equations [36, 37, 38, 39].

Despite the simple structure of the system of equations (1), we can use, for its solution, some elimination methods, symbolic transformations, etc., only in very particular situations. The reason is insufficient algebraic properties of the interval arithmetics. The absence of full distributivity in Kaucher arithmetic makes it generally impossible to perform even such a simple operation as reduction of similar terms. It is for this reason that the methods considered in our work are essentially numerical. An important theoretical result on formal solutions to interval linear systems was obtained by A.V. Lakeyev who managed to show the NP-complexity (intractability) of computing formal (algebraic) solutions to interval systems of linear equations in general form [21, 22].

To find formal solutions of interval linear systems of equations, several numerical methods were proposed, of which the subdifferential Newton method is the most efficient (see [36, 38, 40], and its computer implementations are freely available at [41]). The goal of this paper is the development of stationary single-step iterative methods for computing formal (algebraic) solutions to interval systems of linear algebraic equations. The need to build such methods is due to a number of facts. Despite very high efficacy of the subdifferential Newton method in practice, its justification for the most general case faces a number of difficulties. In addition to calculating the solution, the methods of the type we are going to construct provide also a proof of the uniqueness of the found solution, which is not provided by the subdifferential Newton method. Finally, another reason for the need to develop single-step stationary iterative methods is the fact that they are able to solve interval systems of equations that are not linear in form, for example,

[TABLE]

where $\textbf{{b}}(x)$ is an interval function of the unknown variable $x$.

To date, computational mathematics has accumulated a large arsenal of theoretical approaches and efficient practical algorithms for solving a wide variety of equations and systems of equations. Can we use any of these methods? Is it possible to apply to our problem any of the traditional numerical methods for solving “operator equations”? Yes, but with some reservations and modifications.

The majority of traditional methods for the solution of equations and systems of equations relate to operator equations in linear spaces. Formally, these methods are not applicable to the problem of computing the formal (algebraic) solutions of system (1), since $\mathbb{IR}^{n}$ and $\mathbb{KR}^{n}$ are not linear spaces (see [38]). However, we can easily circumvent this difficulty by embedding the space $\mathbb{IR}^{n}$ or $\mathbb{KR}^{n}$ into the usual well-studied Euclidean space $\mathbb{R}^{2n}$.

As we have already noted, the problem of finding formal solutions to interval equations is, in essence, the traditional mathematical problem of solving some equations, and most of the classical numerical analysis is devoted to solving such problems. But the peculiarity of our situation is that the basic set $\mathbb{KR}^{n}$, on which the equations to be solved are considered, is not a linear space at all: the lack of distributivity in interval arithmetic leads to a violation of the axiom of linear space that requires the fulfillment of the identity

[TABLE]

for all $\textbf{{x}}\in\mathbb{KR}^{n}$ and any scalars $\mu,\nu\in\mathbb{R}$. Thus, most of the existing approaches to the study of operator equations and to the calculation of their solutions are not directly applicable to our problem.

Moreover, remaining within the interval space $\mathbb{KR}^{n}$, we will not be able to perform a theoretical analysis of the situation and understand some phenomena. For example, the point matrix

[TABLE]

is regular (non-singular) in the sense of classical linear algebra, but multiplication by this matrix in $\mathbb{KR}^{2}$ can nullify a non-zero vector:

[TABLE]

What is the reason? It is hardly possible to detect it from within the interval space, which is essentially non-linear. So, there is an urgent need to transfer our considerations to a certain linear space, which we denote by $U$ for generality. We also assume that a topological structure is determined on $U$ consistent with its linear structure.

From an abstract mathematical point of view, we have two different spaces, the interval space $\mathbb{KR}^{n}$ and the linear space $U$, on which essentially different algebraic structures are given. How is it possible to “jump over” from the first one to the second? We are going to do it in a way similar to a change of variables, which is defined in the following subsection and is called immersion.

3.1 Definition and main properties

To transfer our considerations from the interval space to a linear one, we should build some mapping

[TABLE]

— embedding of the interval space $\mathbb{KR}^{n}$ into the linear space $U$. It must be bijective (one-to-one mapping “to”) in order to correctly restore the interval preimage by its image in $U$, and vice versa. Further, it is easy to understand that any bijection $\iota:\mathbb{KR}^{n}\to U$ also generates a bijection from the set of all mappings of $\mathbb{KR}^{n}$ into itself onto the set of all mappings of $U$ to itself. More precisely, each $\varphi:\mathbb{KR}^{n}\to\mathbb{KR}^{n}$ is associated with a uniquely determined mapping

[TABLE]

where “$\circ$” denotes a composition of mappings. Thus, we can argue consider mappings of linear spaces as “exact copies” of interval mappings.

Definition 4* For an interval mapping $\varphi:\mathbb{KR}^{n}\to\mathbb{KR}^{n}$ and a fixed embedding $\;\iota:\mathbb{KR}^{n}\to U$, we will refer to the mapping of the linear space $U$ into itself which is determined by (27) as induced mapping for $\varphi$ (or, expanded, $\iota$-induced).*

Visually, the situation is represented by a commutative diagram in Fig. 1.

The properties of the mappings $\varphi$ and $(\,\iota\circ\varphi\circ\iota^{-1})$ are closely related, so instead of the study of $\varphi$, one can investigate the induced mapping $(\,\iota\circ\varphi\circ\iota^{-1})$. Moreover, we can replace the problem of solving the equation in $\mathbb{KR}^{n}$ with the problem of solving the equation in the linear space $U$, coming to a situation more familiar to modern numerical analysis.

Definition 5* Let the equation be given in the interval space $\mathbb{KR}^{n}$,*

[TABLE]

where $\varphi,\psi:\mathbb{KR}^{n}\to\mathbb{KR}^{n}$ are some mappings, and an embedding $\iota:\mathbb{KR}^{n}\to U$ is fixed. We shall call the induced equation for (28) such an equation

[TABLE]

in the linear space $U$, that $\,\Phi$ and $\Psi$ are induced mappings for $\varphi$ and $\psi$ respectively, i. e. $\Phi=\iota\circ\varphi\circ\iota^{-1}$ and $\Psi=\iota\circ\psi\circ\iota^{-1}$.

Thus, the initial interval equation

[TABLE]

has a formal solution $\textbf{{x}}^{\ast}\in\mathbb{KR}^{n}$ if and only if the induced equation

[TABLE]

has a solution $y^{\ast}\in U$ in the linear space. In this case, the desired formal interval solution $\textbf{{x}}^{\ast}$ for (29) is uniquely reconstructed by $y^{\ast}$ from the relation

[TABLE]

We are interested in the specific situation with interval linear equations (1)–(2). We can change the original problem — finding solution of the equation

[TABLE]

such that

[TABLE]

to the problem of solving equations

[TABLE]

in the linear space $U$ with induced mappings

[TABLE]

defined as

[TABLE]

A more general consideration. Since $\iota$ and $\iota^{-1}$ are bijections, then the invertibility of any mapping $\varphi$ on the interval space is equivalent to the invertibility of the $\iota$-induced map $\Phi:=\iota\circ\varphi\circ\iota^{-1}$ acting on the linear space $U$. Herewith

[TABLE]

The main question concerning the construction of the embedding of the interval space into the linear space is to choose a reasonable compromise between the simplicity of this mapping and the convenient form of induced mappings (27). Among all bijective embedding of $\iota:\mathbb{KR}^{n}\to U$, it makes sense to select special embeddings that

preserve the additive algebraic structure of $\mathbb{KR}^{n}$,

i. e. such that $\iota(\textbf{{u}}+\textbf{{v}})=\iota(\textbf{{u}})+\iota(\textbf{{v}})$ for any $\textbf{{u}},\textbf{{v}}\in\mathbb{KR}^{n}$,

• 2)

preserve the topological structure of $\mathbb{KR}^{n}$,

i. e. such that both the mapping $\iota:\mathbb{KR}^{n}\to U$ itself

and its inverse $\iota^{-1}:U\to\mathbb{KR}^{n}$ are continuous.

Embeddings $\mathbb{KR}^{n}\to U$ satisfying two above prescribed conditions will be called immersions of the interval space $\mathbb{KR}^{n}$ into the linear space $U$. Thus, formally we accept the following

Definition 6* [36]

Let $U$ be a linear space. A bijective mapping $\iota:\mathbb{KR}^{n}\to U$ will be called immersion of $\mathbb{KR}^{n}$ in $U$, if it satisfies the following properties:

   (1) $\iota$ is an isomorphism of additive groups $\mathbb{KR}^{n}$ and $U$,

   (2) $\iota$ is a homeomorphism of topological spaces $\mathbb{KR}^{n}$ and $U$.*

For example, if the interval $\textbf{{v}}\in\mathbb{KR}$ is matched with a pair of numbers $(\underline{\textbf{{v}}},\overline{\textbf{{v}}})\in\mathbb{R}^{2}$, i. e. its endpoints, “forgetting” about their interval sense, then the mapping $\mathbb{KR}\to\mathbb{R}^{2}$ is an immersion. This example is typical in some sense, since, by involving dimension considerations (see, e. g., [7]), it is easy to show that Definition 6 determines the linear space $U$ uniquely: $U$ must be the Euclidean space $\mathbb{R}^{2n}$. This fact is in good agreement with our analytical intuition, and we do not give here its strict substantiation, so as not to overload the overgrown text of the article. The purpose of this preparatory section is the study of the simplest properties of immersions that we will need in the future, when solving the induced equations.

Denote by $0_{\mathbb{KR}^{n}}$ and $0_{\mathbb{R}^{2n}}$ zero vectors in the spaces $\mathbb{KR}^{n}$ and $\mathbb{R}^{2n}$ respectively. It immediately follows from Definition 6, that for any immersion $\iota:\mathbb{KR}^{n}\to\mathbb{R}^{2n}$, we have

[TABLE]

At the same time,

[TABLE]

In addition, the inverse of the immersion map $\iota^{-1}:\mathbb{R}^{2n}\to\mathbb{KR}^{n}$ also satisfies conditions similar to (1)–(2) from Definition 6, and

[TABLE]

Proposition 2*

The immersion is a positive-homogeneous map, i. e.*

[TABLE]

The mapping $\mathbb{R}^{2n}\to\mathbb{KR}^{n}$, inverse to an immersion, is also positively homogeneous.

Proof is standard. Let $\textbf{{x}}\in\mathbb{KR}^{n}$. If $\lambda=k$ is a natural number, then

[TABLE]

If $\lambda=1/l$ for some natural $l$, then, from

[TABLE]

it follows that

[TABLE]

If $\lambda=k/l$ for natural numbers $k$ and $l$, then, using the already considered cases, we obtain

[TABLE]

Hence, the equality $\,\iota(\lambda\textbf{{x}})=\lambda\,\iota(\textbf{{x}})\,$ is valid for all nonnegative rational numbers $\lambda$. Extending it to all non-negative real numbers can be done by passing to the limit, using the continuity of the immersion $\iota$.

The proof for a mapping $\iota^{-1}$ which is inverse to an immersion is performed in a completely similar way. $\blacksquare$

Proposition 3*

If $\iota:\mathbb{KR}^{n}\to\mathbb{R}^{2n}$ is an immersion, and $T$ is a non-singular linear transformation of the space $\mathbb{R}^{2n}$, then $(T\circ\iota)$ is also an immersion.*

Conversely, any other immersion $\varkappa:\mathbb{KR}^{n}\to\mathbb{R}^{2n}$ can be represented as $(T\circ\iota)$ for some nonsingular linear transformation $T:\mathbb{R}^{2n}\rightarrow\mathbb{R}^{2n}$.

Proof. The first part of Proposition is justified trivially.

To prove the second part, we consider the mapping $(\,\varkappa\circ\iota^{-1})$. Being a composition of two isomorphisms, it is an automorphism of the additive group of the linear space $\mathbb{R}^{2n}$, and by virtue of Proposition 3.1 this map is also positively homogeneous. Also, for any $x\in\mathbb{R}^{2n}$

[TABLE]

which implies

[TABLE]

Hence, we get the homogeneity of $(\,\varkappa\circ\iota^{-1})$ with respect to multiplication by negative numbers as well.

Thus, in general, the mapping $\varkappa\circ\iota^{-1}\,$ turns out to be a nonsingular linear transformation of the space $\mathbb{R}^{2n}$. We can therefore take $\,T=\varkappa\circ\iota^{-1}$. $\blacksquare$

3.2 Standard immersion

From Proposition 3.1, it follows that any two immersions of $\mathbb{KR}^{n}$ into $\mathbb{R}^{2n}$ that satisfy Definition 6, are the same to within a linear transformation of $\mathbb{R}^{2n}$. In fact, the choice of a convenient immersion turns out to be even more constrained, since in interval analysis, in addition to the algebraic properties taken into account by Definition 6, the structure of the partial order with respect to inclusion plays an important role. It should also be “adequately” transferred by immersion from $\mathbb{KR}^{n}$ to $\mathbb{R}^{2n}$.

Every immersion $\iota:\mathbb{KR}^{n}\to\mathbb{R}^{2n}$ naturally generates a partial ordering “$\sqsubseteq$” in the linear space $\mathbb{R}^{2n}$, which is the image of the inclusion ordering “$\subseteq$” in $\mathbb{KR}^{n}$ under the immersion $\iota$. Namely,

[TABLE]

Definition 7* The partial ordering “$\,\sqsubseteq$” on $\mathbb{R}^{2n}$, defined through (33), will be called induced partial ordering.*

For any $x,\,y,\,u,\,v\in\mathbb{R}^{2n}$, there holds

[TABLE]

and in such cases it is said that the partial order “$\sqsubseteq$” is consistent with the linear structure on $\mathbb{R}^{2n}$ [12, 42]. As a consequence, this partial order can be equivalently defined by specifying the so-called cone of positive elements, i. e. the set $K_{\sqsubseteq}=\{\,x\in\mathbb{R}^{2n}\mid 0\sqsubseteq x\,\}$ [19, 42].

Recall that a cone in a linear topological space is a closed convex positively invariant set that does not contain any one-dimensional subspace. As is known, in a partially ordered linear space, where the order is consistent with the linear structure, the set of positive elements is a cone. Conversely, the assignment of the cone $K_{\sqsubseteq}$ uniquely determines the partial ordering of the space, for which

[TABLE]

It is clear that the specific formulas determining the induced order of “$\sqsubseteq$” depend on the form of the immersion $\iota$. But on the Euclidean space $\mathbb{R}^{2n}$, the simplest and most convenient is to set the order in a component-wise manner, when

[TABLE]

Accordingly, the cone of positive elements with this ordering $\mathbb{R}^{2n}$ is the set

[TABLE]

i. e., the positive orthant of the space $\mathbb{R}^{2n}$. It is therefore natural to require from the immersion that the order (33) induced by it coincides with this simple component-wise order (34), i. e. that

[TABLE]

For what immersion $\mathbb{KR}^{n}\to\mathbb{R}^{2n}$, is this possible?

It is easy to see that the required immersion is the so-called standard immersion, first introduced in [36]:

Definition 8* The immersion $\,{\mathrm{sti}\;}:\mathbb{KR}^{n}\to\mathbb{R}^{2n}$ that acts according to the rule*

[TABLE]

i. e., such that the negated left endpoints of the intervals $\textbf{{x}}_{1}$, $\textbf{{x}}_{2}$, …, $\textbf{{x}}_{n}$ become the first, second, …, $n$-th component of the point $2n$-vector, while the right endpoints of the intervals $\textbf{{x}}_{1},\textbf{{x}}_{2},\ldots,\textbf{{x}}_{n}$ become the $(n+1)$-st, …, $2n$-th components of the point $2n$-vector respectively will be called standard immersion of the interval space $\mathbb{KR}^{n}$ into $\mathbb{R}^{2n}$.

Corollary. From definition (33) of the induced order on $\mathbb{R}^{2n}$ and the requirement (35) for the standard immersion sti, it is easy to deduce that

[TABLE]

for any bounded family of interval vectors $\{\,\textbf{{x}}_{\vartheta}\in\mathbb{KR}^{n}\mid\vartheta\in\varTheta\,\}$, where $\varTheta$ is an index set. Therefore, the standard immersion maps supremums with respect to inclusion in the interval space $\mathbb{KR}^{n}$ to supremums with the component-wise order on $\mathbb{R}^{2n}$. The same is true for infimums.

So, the fact that the induced partial order on the linear space $\mathbb{R}^{2n}$ coincides with the usual component-wise ordering and, hence, the simplification of calculations and reasoning are the main justification for the form (36) chosen for immersion called “standard”. Moreover, the foregoing rather strongly suggests that in the theoretical part of our work we should consider only the standard immersion of the form (36), although other immersions can sometimes be practically useful. For example, in the computer implementation of the algorithms described in this book, the author often used the simplest immersion

[TABLE]

which is more convenient for practical programming, etc.

It is useful to provide a methodological commentary on the content of this and the preceding subsections. The method of identifying the endpoints of an interval or interval vector with the components of a point vector in the Euclidean space of double dimension was often used by researchers. But we singled out the procedure of this identification into a separate notion — immersion $\mathbb{KR}^{n}\to\mathbb{R}^{2n}$ — and undertook its thorough investigation. For what purpose? Could it be possible to do without “unnecessary abstractions”?

In addition to the fact that explicit and conscious operating with any object is always more preferable than the implicit one, “by default”, there are at least two more reasons to consider the immersion as an independent concept:

• $\bullet$

the mapping $\mathbb{KR}^{n}\to\mathbb{R}^{2n}$ cannot be determined in a unique way once and for all, which would be the most convenient (natural, etc.) for all possible practical situations;

• $\bullet$

we can get a tangible benefit from this non-uniqueness, i. e. the most complete use of the features of immersions in each specific case.

It is easy to see that both of these arguments really apply in our situation.

3.3 Extended multiplier matrix

Theorem 1*

Let the mapping $\,\iota:\mathbb{KR}^{n}\to\mathbb{R}^{2n}$ be an immersion and*

[TABLE]

is operator of multiplication by a point square matrix in the space $\mathbb{KR}^{n}$, i. e.

[TABLE]

for some $Q\in\mathbb{R}^{n\times n}$, $Q=(\,q_{ij})$. Then the induced mapping $\;(\iota\circ\phi\circ\iota^{-1})$ is a linear transformation of the space $\mathbb{R}^{2n}$. For the standard immersion “ sti”, the matrix of this induced linear transformation $(\,{\mathrm{sti}\;}\circ\phi\circ{\mathrm{sti}\;}^{-1})$ is a block $2n\times 2n$-matrix of the form

[TABLE]

where $n\times n$-submatrices $Q^{+}=(\,q_{ij}^{+})$ and $Q^{-}=(\,q_{ij}^{-})$ are positive and negative parts of $Q$, i. e. the matrices made up of positive and negative parts of the elements of $Q$ respectively.

Previously, such matrices were called “concomitant matrices” (see e. g. [39], but then it turned out that this term was already used in mathematics and has a different meaning.

Proof. In order to substantiate the first statement of the Theorem, we need to show the additivity and homogeneity of the mapping

[TABLE]

The additivity of $\phi$ immediately follows from the distributivity relation

[TABLE]

which is valid for any intervals x, $\textbf{{y}}\in\mathbb{KR}$ provided that $q$ is a point (non-interval) value. The immersion $\iota$ and its inverse mapping $\iota^{-1}$ are also additive. Consequently, the induced mapping $\;\iota\circ\phi\circ\iota^{-1}$ is additive as a composition of additive ones.

Next, the operator $\phi$ of multiplication by a matrix is homogeneous, and both the immersion $\iota$ and its inverse mapping $\iota^{-1}$ are positively homogeneous due to Proposition 3.1. Therefore, the composition $\;\iota\circ\phi\circ\iota^{-1}$ is at least positively homogeneous. In addition, for any $x\in\mathbb{R}^{2n}$, there holds

[TABLE]

which proves the homogeneity of the induced mapping $\;\iota\circ\phi\circ\iota^{-1}$ with respect to multiplication on any scalars.

The second statement of the Theorem is a consequence of the definition of the standard immersion (36) and the rule of multiplying a number by the interval

[TABLE]

for which we can give the following convenient equivalent form

[TABLE]

$\blacksquare$

The block $2n\!\times\!2n$-matrix from Theorem 1 is so important in the theory we develop that we will take a special notation and term for it.

Definition 9*

If $Q$ is a point $n\times n$-matrix, we set*

[TABLE]

and will refer to the point $2n\times 2n$-matrix $Q\widetilde{\phantom{p}}$ as extended multiplier matrix for $Q$.

An important feature of the extended multiplier matrix $Q\widetilde{\phantom{p}}\in\mathbb{R}^{2n\times 2n}$ is that it is always non-negative: such a matrix must match a “$\leq$”-isotonic operator on $\mathbb{R}^{2n}$ induced by multiplying by $Q$ in the interval space $\mathbb{KR}^{n}$, which is inclusion isotonic.

Corollary from Theorem 1. Using the definition of the induced map, it is easy to conclude that, for any point $n\times n$-matrix $Q$ and any $x\in\mathbb{R}^{2n}$, there holds

[TABLE]

Similarly, for any point $n\times n$-matrix $Q$ and any interval vector $\textbf{{x}}\in\mathbb{KR}^{n}$, we have

[TABLE]

(The relations (39) and (40) are illustrated by the commutative diagram shown in Fig. 1).

3.4 Absolutely regular matrices

Theorem 2*

For a point $n\times n$-matrix $Q$, the following conditions are equivalent:*

• (a)

for any $x\in\mathbb{KR}^{n}$, $Q\textbf{{x}}=0$  if and only if  $\textbf{{x}}=0$;

• (b)

the matrix $Q\widetilde{\phantom{p}}\in\mathbb{R}^{2n\times 2n}$, extended multiplier matrix for $Q$, is regular;

• (c)

*both the matrix $Q$ and its absolute value $|Q|$ (i. e. the matrix *

formed by moduli of the elements of $Q$) are regular.

Proof. The equivalence “(a) $\Leftrightarrow$ (b)” is a consequence of representation (40).

To prove the equivalence of conditions (b) and (c), we perform the following transformations with the extended multiplier matrix $Q\widetilde{\phantom{p}}$. We add its first row to the $(n+1)$-th row, then we add its second row to the $(n+2)$-th row, and so on, up to the $n$-th row, which we add the $2n$-th row of the matrix $Q\widetilde{\phantom{p}}$. Insofar as

[TABLE]

for any real number $q$, then the result of our transformation will be the following $2n\times 2n$-matrix:

[TABLE]

Next, we subtract the $(n+1)$-th column of the matrix (41) from its first column, then we subtract the $(n+2)$-th column from its second column, and so on, up to the $2n$-th column which we subtract from the $n$-th column of (41). Since

[TABLE]

for any real number $q$, then we get a block-triangular $2n\times 2n$-matrix

[TABLE]

As is known from matrix theory, the transformations we have done (linear combination of rows and columns) do not change the property of a matrix to be regular or singular (see, for example, [8]). Consequently, the matrix (42) is regular or singular along with the extended multiplier matrix $Q\widetilde{\phantom{p}}$. But due to the special form of the matrix (42), its determinant is equal to the product of the determinants of the diagonal blocks, i. e. determinants of the matrices $Q$ and $|Q|$. Thus, the extended multiplier matrix $Q\widetilde{\phantom{p}}$ is regular indeed if and only if both matrices $Q$ and $|Q|$ are regular (non-singular). $\blacksquare$

We have already noted that the regularity (non-singularity) of the point matrix $Q$ in the sense of classical linear algebra does not necessarily imply that the corresponding multiplication operator by $Q$ in $\mathbb{KR}^{n}$ is invertible. But now the phenomenon of the matrix (26) and its like is fully explained: although such matrices themselves may be non-singular, but multiplying by them corresponds, after immersion in linear space, to multiplication by singular extended multiplier matrices. For example, for the matrix (26), the extended multiplier matrix is

[TABLE]

and its determinant is zero!

The result of Theorem 2 allows us to state the following

Definition 10*

A point $n\times n$-matrix $Q$ that satisfies any one (and, hence, every one) of the equivalent conditions (a)–(c) of Theorem 2, is called absolutely regular (absolutely non-singular).*

For example, the identity matrix is absolutely regular, whereas the matrix (26) is regular in the usual sense, but not absolutely regular. It is also obvious that if a matrix is singular in the usual sense, then it is certainly not absolutely regular. All non-negative regular matrices are also absolutely regular. A practically convenient criterion for checking the absolute regularity of a matrix is provided by condition (c) from Theorem 2. For example, instead of calculating the determinant of the extended multiplier matrix (43) for the matrix (26), one could just notice that the moduli matrix

[TABLE]

is singular.

The class of matrices, determined by the condition of Definition 10 and equivalent to it, was introduced by the author in the works of the past century [36, 38]. But the term “absolutely regular matrix” was coined relatively recently, already in the 2000s. Previously, these matrices were called $\iota$-regular ($\iota$-nonsingular), and then completely regular (completely nonsingular) matrices [39]. A long and nontrivial way indeed, but the term “absolutely regular matrix” seems to be the most adequate to the meaning of the concept.

Corollary from Theorems 1 and 2. The operator $\phi:\mathbb{KR}^{n}\to\mathbb{KR}^{n}$ determined by the multiplication by a square point matrix in $\mathbb{KR}^{n}$, i. e. such that

[TABLE]

is invertible if and only if the matrix $Q$ is absolutely regular. Then the inverse operator $\phi^{-1}:\mathbb{KR}^{n}\to\mathbb{KR}^{n}$ acts, according to (30), as follows:

[TABLE]

Despite the existence of the explicit formula (44), an operator that is inverse to the operator of multiplication by a point $n\times n$-matrix $Q$ in $\mathbb{KR}^{n}$ cannot generally be expressed through multiplication by a matrix from $\mathbb{KR}^{n}$ (in particular, by the matrix $Q^{-1}$). This follows, for example, from the fact that the matrix $(Q\widetilde{\phantom{p}})^{-1}$ does not need to be non-negative, and then the inverse operator is not monotonic with respect to inclusion in $\mathbb{KR}^{n}$. But multiplication by a point matrix in $\mathbb{KR}^{n}$ is always monotone with respect to inclusion.

The above Corollary is formulated in an abstract manner, but it has quite practical applications. Using formula (44), one can easily find formal solutions to interval systems of linear equations in which the square matrix is a point absolutely regular matrix, and only the right-hand side vector is interval, i. e. of the form

[TABLE]

Consider, for example, the interval linear system

[TABLE]

First of all, we note that the matrix of this system is point (non-interval) and absolutely regular, since both the matrix itself and its modulus are regular:

[TABLE]

Then the multiplication operator by this matrix is invertible in the interval space $\mathbb{KR}^{n}$, and its inverse is given by formula (44). To apply it for computing formal solutions to (46), we construct the extended multiplier matrix

[TABLE]

and then multiply its inverse by the result of the standard immersion “sti” acting on the right-hand side vector of (46), that is, by $(0,-10,10,20)^{\top}$:

[TABLE]

Then the formal solution for system (46) is

[TABLE]

By direct substitution, we can verify that the resulting interval vector is indeed a formal solution of system (46).

S. Markov proposed in [25] an equivalent method for finding formal solutions to interval linear systems with point (non-interval) matrices, i. e. of the form (45). Its essence, if we give up unnecessary abstractions, is as follows.

If $\textbf{{x}}^{\ast}\in\mathbb{KR}^{n}$ is a formal solution to system (45), then, taking the midpoint and radius of both sides of the equality

[TABLE]

we get

[TABLE]

Thus, computing the formal solution is reduced to the solution of two ordinary systems of linear equations with the matrices $A$ and $|A|$, from which we find the middle of the formal solution $\textbf{{x}}^{\ast}$ and its radius. It is clear that a sufficient condition for the existence of these solutions is the non-singularity of both matrices $A$ and $|A|$. If we want the solvability of these systems of equations for the middle and radius of the formal solution for any right-hand sides b, then the matrix $A$ must be absolutely regular.

4 The subdifferential Newton method

In this section of the paper, we briefly present results on the subdifferential Newton method for finding formal solutions of interval linear systems.

Recall some definitions and facts from convex analysis [24, 33]. Let $U,V$ be real linear spaces, and $V$ be an ordered linear space with an order “$\preceq$”. The mapping $F:U\rightarrow V$ is called order convex with respect to “$\preceq$”, if

[TABLE]

for any $x,\,y\in U$ and $\lambda\in\;]0,1[$. Subdifferential of the order convex mapping $F:U\rightarrow V$ at the point $x$ is the set $\partial_{\preceq}F(x)$ of all linear operators $D:U\rightarrow V$ such that

[TABLE]

for any $y\in U$. It is known that for each interior point $x$ from the effective domain of definition of the order convex function $F$, the subdifferential is a non-empty set, whose elements are called subgradient of $F$ at $x$.

Suppose that the interval matrix A is such that in each of its rows all elements are either only proper or only improper, and let the sets of natural numbers $I^{\prime}=\{\,i^{\prime}_{1},i^{\prime}_{2},\ldots,i^{\prime}_{\alpha}\,\}$ and $I^{\prime\prime}=\{\,i^{\prime\prime}_{1},i^{\prime\prime}_{2},\ldots,i^{\prime\prime}_{\beta}\,\}$, $I^{\prime}\cap I^{\prime\prime}=\varnothing$, $\alpha+\beta=n$, represent row indices of A with proper and improper intervals, respectively:

[TABLE]

If we define on $\mathbb{R}^{2n}$ a partial order “$\sqsubseteq$” as follows

[TABLE]

then, as a consequence of the sub- and superdistributive properties of Kaucher arithmetic, the induced mapping $\Phi(x)$, determined by (27) for the interval mapping $\textbf{{x}}\mapsto\textbf{{Ax}}\ominus\textbf{{b}}$, is order convex with respect to “$\sqsubseteq$”. Accordingly, at any point $x\in\mathbb{R}^{2n}$ there will be a non-empty subdifferential $\partial_{\sqsubseteq}\Phi(x)$, easily computable, because $\Phi$ is a polyhedral (piecewise affine) mapping. The subdifferential Newton method proved to be an efficient tool for finding formal solutions to such interval linear algebraic systems. It is a further development of the well-known results on monotonically convergent Newton-type methods in ordered linear spaces.

The subdifferential Newton method with a special starting vector

 

As a starting approximation of $x^{(0)}\in\mathbb{R}^{2n}$, we take the solution of the “midpoint” $2n\!\times\!2n$-system

$(\mathrm{mid}\;\textbf{{A}})\widetilde{\phantom{p}}\,x={\mathrm{sti}\;}(\textbf{{b}}).$

If the $k$-th approximation $x^{(k)}\in\mathbb{R}^{2n}$, $k=0,1,\ldots\,$, is already known, then we find some subgradient $D^{(k)}\in\partial_{\sqsubseteq}\Phi(\,x^{(k)})\,$, $D\in\mathbb{R}^{2n\times 2n}$, and then take

$x^{(k+1)}\leftarrow x^{(k)}-\tau\left(\,D^{(k)}\right)^{-1}\left(\Phi(\,x^{(k)})\right).\vskip 12.0pt plus 4.0pt minus 4.0pt$

In the above pseudocode, $\tau$ is a damping factor from $]0,1]$. There holds

Theorem 3* If, in each row of the interval matrix A, all elements are only proper or only improper, the proper projection ${\mathrm{pro}\;}\textbf{{A}}$ contains only absolutely regular point matrices and is sufficiently narrow (i. e. $\|\mathrm{rad}\;\textbf{{A}}\|$ is small enough), then the subdifferential Newton method converges in $\mathbb{R}^{2n}$ to ${\mathrm{sti}\;}(\textbf{{x}}^{*})$, where $\textbf{{x}}^{*}$ is a formal solution of the system $\textbf{{Ax}}=\textbf{{b}}$ and “sti” is the standard immersion of $\mathbb{KR}^{n}$ to $\mathbb{R}^{2n}$.*

As applied to more particular types of interval systems, the proof of the convergence of the subdifferential Newton method and its detailed discussion are contained in [38]. In practice, the subdifferential method of Newton converges extremely quickly, in a small finite number of iterations (especially for $\tau$ close or equal to one), which is explained by the polyhedral nature of the induced function $\Phi(x)$ that corresponds to interval linear systems. Another advantage of this method is the absence of any problems with the choice of the starting vector. Finally, the subdifferential Newton method works well in practice even for such interval systems in which proper and improper elements are arbitrarily mixed in the interval matrix. In general, the Newton subdifferential method is of great practical importance, but in this paper we will consider it mainly as an auxiliary algorithm, the initial link in the composition of numerical processes with a wider domain of convergence.

5 Stationary single-step iterative methods

5.1 General approach: distributive splittings

How else can numerical methods be proposed for finding formal (algebraic) solutions to interval linear systems of equations (1)–(2) with arbitrary matrices and right-hand sides?

Although interval arithmetic differs from the real axis $\mathbb{R}$ on its algebraic properties, the simple structure of system (1)–(2) makes it possible to apply for this purpose some general ideas related to stationary iterative methods for ordinary linear systems of equations. In particular, we can construct interval versions of Jacobi method, Gauss-Seidel method, etc. When implementing such methods, it is sometimes possible to iterate directly in the interval space $\mathbb{IR}^{n}$ or $\mathbb{KR}^{n}$, even without immersing it in the linear space of double dimension $\mathbb{R}^{2n}$. On the other hand, the convergence obtained in such stationary iteration algorithms is principally slower than that for the subdifferential (quasidifferential) Newton method.

It should be noted that iterative Jacobi-type methods for finding formal solutions to interval linear systems of equations have already been considered by different authors. V. Zyuzin was a pioneer in this area, who proposed them in the late 80s of the last century [45, 46, 47]. Then these methods were used in the works of L. Kupriyanova [20], S. Markov [25, 26], and Spanish researchers [34, 35]. But the applicability scope and convergence rate of these simplest methods are unsatisfactory, which necessitates the development of more advanced numerical methods. Below, we consider a more general construction than Jacobi-type methods.

In accordance with the general scheme of stationary single-step iterative methods, the original equation (1)–(2) must be equivalently reduced to a fixed-point form

[TABLE]

with some operator $\textbf{{T}}:\mathbb{KR}^{n}\rightarrow\mathbb{KR}^{n}$. Then, after selecting an initial approximation $\textbf{{x}}^{(0)}$, iterations are started:

[TABLE]

Under some special conditions on the transition operator T (when it is a contraction, etc.) and on the initial approximation $\textbf{{x}}^{(0)}$, the sequence $\textbf{{x}}^{(k)}$ converges to a fixed point of T, i. e. to the desired formal solution of the equation (2). But, in contrast to the traditional non-interval case, bringing the interval system (1)–(2) to the form (47) is not a trivial task due to insufficient algebraic properties of Kaucher arithmetic. The problem is that at least two terms with the variable x in the formula (47) (which is equivalent to $\,\textbf{{x}}\ominus\textbf{{T}}(\textbf{{x}})=0\,$) should eventually collapse into the expression $\,\textbf{{Ax}}\ominus\textbf{{b}}$, containing only one occurrence of the variable x. In the absence of a full-fledged opportunity to reduce such terms, special means are required for transforming the original equation (2) into the form (47).

One of the promising approaches to the construction of iterative schemes for the solution of the equation (2) is to go the “opposite way”, from possible representations

[TABLE]

where

• (i)

the function $\mathcal{G}:\mathbb{KR}^{n}\to\mathbb{KR}^{n}$ is “easily invertible”, i. e. its inverse function $\mathcal{G}^{-1}:\mathbb{KR}^{n}\to\mathbb{KR}^{n}$, such that $\mathcal{G}^{-1}(\,\mathcal{G}(\textbf{{x}})\,)=\ =\mathcal{G}(\,\mathcal{G}^{-1}(\textbf{{x}})\,)=\textbf{{x}}$, can be easily constructed for $\mathcal{G}$;

• (ii)

the function $\mathcal{H}:\mathbb{KR}^{n}\to\mathbb{KR}^{n}$ is easily computable.

Definition 11*

Let $\psi:\mathbb{KR}^{n}\to\mathbb{KR}^{n}$ be an operator in $\mathbb{KR}^{n}$ determined by multiplication by the interval matrix A, i. e. such that $\,\psi(\textbf{{x}})=\textbf{{Ax}}$. The representation of $\psi$ in the form of (49), satisfying for any $\textbf{{x}}\in\mathbb{KR}^{n}$ the above conditions (i)–(ii) will be called a splitting of the operator $\psi$ of multiplication by the matrix A, or, briefly, splitting of the matrix A.*

If we know some splitting of matrix A in the interval linear system

[TABLE]

then we can proceed to the equivalent equation

[TABLE]

or

[TABLE]

which coincides with the desired form (47). Accordingly, the iterative process for computing the formal solution can be organized by the formula

[TABLE]

after chosing a starting approximation $\textbf{{x}}^{(0)}$. Further, we restrict ourselves to the simplest case where the functions $\mathcal{G}$, $\mathcal{H}:\mathbb{KR}^{n}\to\mathbb{KR}^{n}$ themselves either represent multiplications by some interval matrices, or they are mappings similar in form.

Suppose that, in system (2), the interval matrix A is regular. Then there are at least two possibilities for the splitting (49), i. e. for the splitting of the multiplication operator by A, with an easily invertible mapping $\mathcal{G}^{-1}(\cdot)$:

• A)

$\mathcal{G}(\cdot)$ is taken in the form of multiplication by a point absolutely regular matrix $G$, that is,

[TABLE]

In this case, $\,\mathcal{H}(\textbf{{x}})=\textbf{{Hx}}$, $\textbf{{H}}=\textbf{{A}}-G$, and nonzero elements of $G$ and H are chosen so that they have the same signs (due to this fact, equality (49) is ensured for all $\textbf{{x}},\textbf{{y}}\in\mathbb{KR}^{n}$ by virtue of distributivity (24)). The inverse mapping $\mathcal{G}^{-1}(\cdot)$ is then determined according to Corollary from Section 3.4.

• B)

$\mathcal{G}(\cdot)$ and $\mathcal{H}(\cdot)$ are taken, respectively, in the form

[TABLE]

where G and H are upper and lower triangular interval matrices, $\textbf{{A}}=\textbf{{G}}+\textbf{{H}}$, where G has invertible elements on the main diagonal, and H has zero main diagonal (perhaps, to do this, one first need to swap the equations of the system).

Then the inverse mapping of $\mathcal{G}^{-1}(\cdot)$ is such that the result y of its action on the element $\textbf{{x}}\in\mathbb{KR}^{n}$ is determined by the “backward substitution” formulas for the triangular system Gy = x. It is natural to use the term triangular splitting with respect to such a splitting of the multiplication operator by A.

Note that in both cases the inverse mapping $\mathcal{G}^{-1}:\mathbb{KR}^{n}\to\mathbb{KR}^{n}\,$, generally speaking, cannot be defined through multiplication by an interval matrix.

5.2 Splitting an absolutely regular matrix

In this section, for the first of the matrix splitting cases considered in §5.1, we write out the computation formulas of the corresponding iterative process (50) in an explicit form. We also show how it is possible in practice to construct a splitting of an interval matrix.

Definition 12* For $\textbf{{a}}\in\mathbb{KR}$, we denote*

[TABLE]

that is, taking the element which is closest to zero from the proper projection of the interval a.

$\lfloor\textbf{{a}}\rfloor$ is the point with the smallest absolute value from the proper projection of the interval, and it has the same sign as the interval itself. It is easy to understand that if a real number $a$ lies between [math] and $\lfloor\textbf{{a}}\rfloor$, i. e. $a\in 0\lor\lfloor\textbf{{a}}\rfloor$, the intervals $(\textbf{{a}}-a)$ and $[a,a]$ have the same sign, and therefore constitute a splitting for the multiplication operator on a in the sense of Definition 11. Therefore, the condition A of §5.1 can be satisfied, for example, if we take

[TABLE]

where $\lfloor\textbf{{A}}\rfloor$ means element-wise application to A of the unary operation $\lfloor\cdot\rfloor$, and the interval matrix in the right-hand side of the inclusion is obtained by using the operation (9).

To minimize the absolute value of the remainder $\textbf{{H}}=\textbf{{A}}-G$, one can assign $G=\lfloor\textbf{{A}}\rfloor$. It is clear that, under the assumption we made about the regularity of A, the matrix $G$ also turns out to be non-singular. If, additionally, $G$ is absolutely regular, then the inverse mapping $\mathcal{G}^{-1}(\cdot)$ corresponds to multiplication by the matrix $(G\widetilde{\phantom{p}})^{-1}$ in $\mathbb{R}^{2n}$ (see Section 3.4). In any case, we always have the opportunity to make the matrix $G$ absolutely regular by slightly reducing the absolute value of its nonzero elements that do not violate the splitting condition (51).

We write out the formulas of the resulting iteration process in the Euclidean space $\mathbb{R}^{2n}$. As a specification for (50), we get

[TABLE]

where “sti” is the standard immersion of $\mathbb{KR}^{n}$ into $\mathbb{R}^{2n}$ and

[TABLE]

The iterative process with such splitting works satisfactorily, but sometimes not as well as one would like. For example, it does not lead to success in solving the Barth-Nuding interval linear system [1],

[TABLE]

Therefore, it makes sense to consider other recipes for splitting the matrix of the interval linear system.

Yet another way of splitting the interval matrix can be based on the generalized distributivity law (25) suggested by S. Markov. We introduce

Definition 13* For $\textbf{{a}}\in\mathbb{KR}$, we denote*

[TABLE]

that is, taking the largest, in absolute value, element from the proper projection of the interval, if it does not contain zero, and zero otherwise.

If $0\not\in{\mathrm{pro}\;}\textbf{{a}}$, then $\lceil\textbf{{a}}\rceil$ is the point from the proper projection of the interval a, which has the largest absolute value and the same sign as the interval itself (unlike $\lfloor\textbf{{a}}\rfloor$). It is easy to see that if a real number $a$ coincides in sign with $\lceil\textbf{{a}}\rceil$ and $|a|>\lceil\textbf{{a}}\rceil$, then the intervals $(\textbf{{a}}-a)$ and $[a,a]$ have different signs, and the sign of their sum $(\textbf{{a}}-a)+a$ is the same as for a. Therefore, for any $\textbf{{x}}\in\mathbb{KR}$, by virtue of (25), there holds

[TABLE]

Consequently, in the general formula of the iterative processes (50), we can put

[TABLE]

In $\mathbb{R}^{2n}$, an explicit formula of the iterative process based on the above point splitting of the matrix of the interval system has the form

[TABLE]

where the matrix $G\in\mathbb{R}^{n\times n}$ and the mapping $\mathcal{H}(\cdot)$ are defined by (55)–(57), and “sti” is the standard immersion of $\mathbb{KR}^{n}$ in $\mathbb{R}^{2n}$. Below, in Section 6, we present the results of numerical experiments with this method appearing under the name ARMSplit (derived from Absolutely Regular Matrix Splitting), which show that it works significantly better than the process (52).

What are the conditions for convergence of the considered iteration processes? The following result is valid:

Theorem 4*

Let the matrices $G\in\mathbb{R}^{n\times n}$ and $\textbf{{H}}\in\mathbb{KR}^{n\times n}$ be the result of splitting, either (53) or (55)–(57), applied to the interval matrix A, and $\mathcal{V}$ is the $2n\times 2n$-matrix $(G\widetilde{\phantom{p}})^{-1}$. If the spectral radius of the matrix $|\mathcal{V}|\,|\textbf{{H}}|\widetilde{\phantom{p}}$ is less than one,*

[TABLE]

then a formal solution $\textbf{{x}}^{\ast}$ of the interval linear system $\textbf{{A}}x=\textbf{{b}}$ exists and is unique, and iterations (52) and (58) converge in $\mathbb{R}^{2n}$ to ${\mathrm{sti}\;}(\textbf{{x}}^{*})$, where “ sti” is the standard immersion of $\mathbb{KR}^{n}$ to $\mathbb{R}^{2n}$.

Before proceeding to the proof of the theorem, let us recall some facts and concepts. In order to consider convergence in the interval spaces and in usual Euclidean spaces, we need to have a “distance”, a “measure of the deviation” of one vector from another. For ordinary linear spaces, distance can be introduced by various equivalent ways that are well known. The most popular is the concept of a metric (distance), as a non-negative real-valued function, which associates with each pair of vectors a non-negative number, the distance between them.

The distance between intervals in the one-dimensional case is known to be defined as

[TABLE]

In the same way, one can determine the distance in the multidimensional case, for interval vectors, but another construction is more popular.

In multidimensional spaces, it is often convenient to work with a vector-valued metric (a vector-valued distance), often called multimetric. For spaces of interval vectors $\mathbb{IR}^{n}$ and $\mathbb{KR}^{n}$, it is natural to introduce it as

[TABLE]

that is, as a vector of distances between the components of the vectors a and b. The multimetric is also defined in a similar way for ordinary linear vector spaces. It is easy to show that it has all the properties of a traditional scalar distance, i. e. non-negativity, symmetry, triangle inequality. With respect to multimetrics, we can talk about convergence, topological completeness and other important and necessary concepts.333Multimetric is called “pseudometric” in the popular book by L. Collatz [4]. We do not follow this usage because, in modern mathematics, the term “pseudometric” has a different meaning. It is “a distance on a set that fulfils the same properties as a metric except relaxes the definition to allow the distance between two different points to be zero” (a citation from MathWorld – A Wolfram Web Resource).

Further, we will substantially rely on the following property of interval multimetrics. For any interval matrix $\textbf{{A}}=(\,\textbf{{a}}_{ij})\in\mathbb{KR}^{m\times n}$ and any interval vectors $\textbf{{x}}=(\,\textbf{{x}}_{j})$, $\textbf{{y}}=(\textbf{{y}}_{j})\in\mathbb{KR}^{n}$, there holds

[TABLE]

Indeed, by virtue of the algebraic properties of interval operations, we can conclude that

[TABLE]

for all $i=1,2,\ldots,m$, which proves the multidimensional estimate (61).

Proof. The convergence of the iterative processes (52) and (58) in the Euclidean space $\mathbb{R}^{2n}$ will be proved in some specially designed multimetric. Namely, we introduce a multimetric $\zeta:\mathbb{R}^{2n}\to\mathbb{R}^{2n}_{+}$ on $\mathbb{R}^{2n}$ as follows:

[TABLE]

Recalling definition (36) of the standard immersion “sti”, we can determine the multi-metric $\zeta$ in another way, specifically

[TABLE]

The proof of Theorem 4 will be first conducted for the iterative process (52). We show that, with respect to the above defined multi-metric $\zeta$, the transition operator $T$ of the iteration scheme (52), determined as

[TABLE]

satisfies the conditions of the Schröder fixed-point theorem (see e. g. [4, 30]).

There holds

[TABLE]

Using inequality (61), we can continue our calculations as follows:

[TABLE]

Therefore,

[TABLE]

From the Frobenius formulas for inverting block matrices (see, for example, [8]), it follows that the $2n\times 2n$-matrix $\mathcal{V}$ is a block matrix of the same structure as the extended multiplier matrix $G\widetilde{\phantom{p}}$, i. e., it is divided into four $n\times n$ blocks, and the diagonal blocks are the same. Therefore, in general, we obtain

[TABLE]

as required.

It is easy to see that the proof conducted is easily adaptable for the iterative process (58) as well, since $|\mathcal{H}(\textbf{{x}})|=|\textbf{{Hx}}|$ for any $\textbf{{x}}\in\mathbb{KR}^{n}$. $\blacksquare$

The main concern of the developers of iterative methods of the form (48) is to reduce, as much as possible, the spectral radius (or norm) of the so-called Lipschitz operator for the transition operator T, in order, first, to ensure the convergence of iterations, and secondly, to accelerate this convergence where it already exists. As follows from the proof of Theorem 4, the matrix of this Lipschitz operator is $|\mathcal{V}|\,|\textbf{{H}}|\widetilde{\phantom{p}}$ for the iterations (52). Optimization of the distributive splitting of the matrix A into $G$ and H is an interesting, but difficult task, and here we will not consider its solution in the most general form. We only note that splitting a point addend is especially convenient in practice for cases where the matrix of the interval systems has many point (non-interval) elements, while the proportion of substantially interval elements in the matrix is small.

5.3 Triangular splitting of the matrix of the equation system

With a triangular splitting of the interval matrix A, the identity

[TABLE]

is achieved due to the fact that the matrices G and H form a disjoint decomposition for A, since nonzero elements in G and H are mutually exclusive. The pseudocode of the iteration process (50) in $\mathbb{KR}^{n}$, with the triangular splitting of the matrix A, which we call TrnSplit, has the form presented in Table 2, where “$\oslash$” is an internal division in $\mathbb{KR}$, i. e. multiplication by the inverse interval (if it exists).

It is worth noting that the iterative method with the triangular splitting of the form we used, is called the Gauss-Seidel method for ordinary linear systems (see e. g. [43]). It would be logical to call our algorithm the “interval Gauss-Seidel method”, but this name has already been fixed with respect to another iterative process (see e. g. [28, 40]). Therefore, we call our method simply a “method based on triangular splitting”.

The study of the convergence of the algorithm TrnSplit was carried out by A.Yu. Karlyuk under the guidance of the author in [13], and the main result of this article was as follows:

Theorem 5* For the interval matrix A of the equation system (2), let the point $n\times n$-matrices $D$, $L$, $R$ are defined by the formulas*

[TABLE]

where “inv” means taking the inverse element with respect to multiplication in Kaucher interval arithmetic. If the matrix

[TABLE]

satisfies $\rho(\mathcal{Q})<1$, then the iteration process TrnSplit converges to a unique formal solution $\textbf{{x}}^{\ast}$ of system (2) from any starting approximation $\textbf{{x}}^{(0)}$. Additionally, the following estimate is valid

[TABLE]

where the multimetrics “ Dist” is defined in (60).

For $\rho\,(\mathcal{Q})<1$ to take place, it suffices, for example, to satisfy the following condition on the interval matrix A of system (2): recurrently calculated numbers $s_{1}$, $s_{2}$, …, $s_{n}$, such that

[TABLE]

are all strictly less than one. In turn, conditions (64) are certainly fulfilled for interval matrices with the property of diagonal dominance:

[TABLE]

A classical result of R. Varga [43] states that, in the point case of $\textbf{{A}}=G-H=C$, $\textbf{{b}}=d$, the method (50) converges the faster the smaller $H$ is.

6 Numerical examples

**Example 1 ** Let us consider the interval linear system by Barth-Nuding [1]

[TABLE]

The algorithm ARMSplit after 10 iterations gives 3 valid significant digits of the exact answer $(-\frac{1}{3},\frac{1}{3})^{\top}$, and after 20 iterations it gives 6 correct significant digits, which in order of labor costs is comparable to the iterative method from [45, 46, 47], based on the cleavage of the main diagonal of the matrix of the interval system. The same figures are achieved by the algorithm ARMSplit when finding a formal solution of system (54) with the dualized matrix (which occurs, for example, in the internal evaluation of the united solution set for (54)). $\blacksquare$

**Example 2 ** Consider an interval linear $40\times 40$-system with the matrix

[TABLE]

and right-hand side vector

[TABLE]

The matrix (65) is obtained from a popular tridiagonal matrix approximating the second derivative on a uniform grid by 10% broadening of the elements, and the right-hand side of the system is obtained by the same broadening of the vector $(\,1,2,3,\ldots,39,40\,)^{\top}$. As in the case of the subdifferential Newton method, neither this system of equations nor the system derived from it by matrix dualization represent a serious problem for the algorithms presented in this paper. The method ARMSplit based on point splitting finds 12–13 correct significant digits for the endpoints of the components of formal solutions of both the original interval system and the system with a dualized matrix after 16 iterations. $\blacksquare$

**Example 3 ** Consider an interval linear $40\times 40$-system with Neumaier matrix [28] having the number $40$ at the main diagonal, i. e.

[TABLE]

and the right-hand side vector

[TABLE]

The algorithm ARMSplit computes, after 40 iterations, the approximation to the formal solution

[TABLE]

with an accuracy of about $10^{-8}$.

The same can be observed when calculating the formal solution to the interval linear system with the matrix dualized to (67) and the right-hand side (68). Note that in this example, the interval matrix of the system is singular (see [28]), but this does not prevent the algorithm ARMSplit from successfully finding a formal solution. $\blacksquare$

**Example 4 ** For the interval linear $7\times 7$-system

[TABLE]

from the work [36], the algorithm ARMSplit diverges, but its formal solution can be successfully found with the use of the subdifferential Newton method (after 9 iterations and for the damping factor $\tau=1$).

When the (7,7)-element of the interval matrix is narrowed, the convergence of the algorithm ARMSplit to the formal solution appears, but it is very slow. For example, with $\textbf{{a}}_{77}=[8,82]$, the algorithm requires about a hundred iterations to get 5 correct significant digits of the answer. $\blacksquare$

Summarizing the last example, we can say that it demonstrates the advantage of the subdifferential Newton method over stationary iterative methods not only in terms of efficiency, but also in terms of the applicability scope.