Computing best discrete least-squares approximations by first-degree splines with free knots

TL;DR

This paper introduces an algorithm for finding the optimal least-squares approximation of discrete functions using first-degree splines with free knots, ensuring a global best fit after finite steps, with practical applications demonstrated.

Contribution

The paper presents a novel finite-step algorithm for computing globally optimal first-degree spline approximations with free knots for discrete data.

Findings

Algorithm guarantees global optimality after finite steps

Numerical examples demonstrate practical effectiveness

Applications include medical data approximation

Abstract

We present an algorithm to compute best least-squares approximations of discrete real-valued functions by first-degree splines (broken lines) with free knots. We demonstrate that the algorithm delivers after a finite number of steps a (global) best approximation. The analysis is complemented by remarks on programming and by a number of numerical examples including applications from medicine (MBC, MIC).

Tables

Table 1. Table 1: Number of regular position vectors as a function of knots and data points.

number of	number of data points $\mu +2$
knots $k$	15	20	25	30	35	40
1	25	35	45	55	65	75
2	265	545	925	1405	1985	2665
3	1561	4991	11521	22151	37881	59711
4	5641	29961	97281	241601	506921	947241
5	13073	124515	590557	1937199	5060441	11326283
6	19825	369305	2668525	11847485	39146185	106114625
7	19825	795455	9173505	56610575	249673265	799538175

Equations

$$S^m(t_1,\dots ,t_k):=\left\{s\in C^{m-1}[a,b]|s_j{|}_{(t_j,t_{j+1})}\in {\mathbb{P}}_m,j=0,1,\dots ,k\right\}$$

$$S_k^m[a,b]$$

$${\Vert F-s^{\ast }\Vert }_{2,X}:={\Vert F-s^{\ast }(X)\Vert }_2=\underset{s\in S_k^1[a,b]}{\inf }{\Vert F-s(X)\Vert }_2.$$

$$0\le p(j)\le p(j+1)\le 2\mu ,j=1,2,\dots ,k-1,$$

$$p(j)=p(j+1)\Rightarrow \text{p(j) is even},$$

$$p(j)\;:=\;\left\{\begin{array}[]{lcl}2i-1&,&\text{if}\;\;t_{j}=x_{i}\\ 2i&,&\text{if}\;\;t_{j}\in\left(x_{i},x_{i+1}\right)\quad.\end{array}\right.$$

$$p(j)=p(j+1)\iff x_{\frac{p(j)}{2}}<t_j<t_{j+1}<x_{\frac{p(j)}{2}+1}.$$

$$x_{\frac{p(j)}{2}}<t_j<t_{j+1}<\dots <t_{j+q}<x_{\frac{p(j)}{2}+1}$$

$$t_j:=x_{\frac{p(j)+1}{2}}\text{for j=1,2,…,k}.$$

$$i_l:=\frac{p\left(j_l\right)}{2},l=1,2,\dots ,r.$$

$$S_l:=\left\{x_{i_l+1},x_{i_l+2},\dots ,x_{i_{l+1}}\right\}.$$

$$t_{j_l+i}=x_{\frac{p\left(j_l+i\right)+1}{2}},i=1,2,\dots ,q_l,$$

$$t_i=x_{\frac{p(i)+1}{2}},i=1,2,\dots ,j_1-1.$$

$$t_i=x_{\frac{p(i)+1}{2}},i=j_r+1,\dots ,k.$$

$$\frac{p(j+l+1)}{2}-\frac{p(j)}{2}\ge 2+l.$$

$$\frac{p(j+1)}{2}-\left(\frac{p(j)}{2}+1\right)+1\ge 2\iff p(j+1)-p(j)\ge 4.$$

$$\frac{p(j)}{2}+1<\frac{p(j+1)+1}{2}\iff p(j+1)-p(j)>1\iff p(j+1)-p(j)\ge 2.$$

$$p(j+1)-p(j)\ge 2$$

$$\big{\{}p(1),\ldots,p(k)\big{\}}\subseteq\big{\{}1,2,\ldots,2\mu-1\big{\}}\quad\text{and}\quad p(1)<p(2)<\ldots<p(k)\quad.$$

$$p(j)\in\big{\{}1,3,5,\ldots,2\mu-3,2\mu-1\big{\}}\;,\;\;j=1,\ldots,k\quad.$$

$$s_{l}\in S^{1}\Big{(}t_{j_{l}+1},\ldots,t_{j_{l}+q_{l}}\Big{)}\;\;=\;\;S^{1}\left(x_{\frac{p\left(j_{l}+1\right)+1}{2}},\ldots,x_{\frac{p\left(j_{l+1}-1\right)+1}{2}}\right)$$

$$s(x):=\left\{\begin{array}[]{lcl}s_{0}(x)&,&x_{0}\leq x\leq z_{1}\\ s_{l}(x)&,&z_{l}<x\leq z_{l+1}\;,\;1\leq l\leq r-1\\ s_{r}(x)&,&z_{r}<x\leq x_{\mu+1}\hskip 85.35826pt.\end{array}\right.$$

$$\big{\|}f-s^{\ast}\big{\|}_{2,X}\;\;=\;\;\inf_{s\in S^{1}_{k}[x_{0},x_{\mu+1}]}\big{\|}f-s\big{\|}_{2,X}\quad.$$

$${\Delta }_j^{(1)}:=\frac{f_{j+1}-f_j}{x_{j+1}-x_j},0\le j\le \mu$$

$${\Delta }_j^{(2)}:=\frac{f_{j+1}-2f_j+f_{j-1}}{{\left(x_{j+1}-x_j\right)}^2},1\le j\le \mu .$$

$$\begin{array}[]{l|lllll|r}&t_{1}&t_{2}&t_{3}&t_{4}&t_{5}&\left\|F-s_{k}(X)\right\|_{2}\\ \hline\cr\hline\cr k=1&8&&&&&0.87586\\ k=2&7&8&&&&0.78881\\ k=3&7&8&9&&&0\\ k=4&4&7&8&9&&0\\ k=5&1&2&7&8&9&0\end{array}$$

$$\begin{array}[]{r||r||c|rrrrr}k&\|F-s_{k}(X)\|_{2}&i&1&2&3&4&5\\ \hline\cr\hline\cr 3&0.2632&t_{i}&858.4883&897.8327&940.2917&&\\ &&s_{k}(t_{i})&0.7642&2.3065&0.6659&&\\ \hline\cr 4&0.1875&t_{i}&831.4392&866.8552&897.5429&940.2917&\\ &&s_{k}(t_{i})&0.7074&1.0396&2.3177&0.6659&\\ \hline\cr 5&0.1349&t_{i}&831.4392&866.8552&898.3019&930.6129&958.3397\\ &&s_{k}(t_{i})&0.7074&1.0396&2.3494&0.9535&0.6153\\ \end{array}$$

$${\phi }_j:=\frac{\text{number of surviving germs after exposition}}{\text{total number of germs before exposition}}$$

$$\left(0,{\phi }_0\right),\left(1,{\phi }_1\right),\dots ,\left(z,{\phi }_z\right).$$

$$\text{MBC}:=\kappa \left(t_1\right)={\kappa }_0\cdot 2^{-t_1}$$

Full text

Computing best discrete least-squares approximations by first-degree splines with free knots

Ludwig J. Cromme111Numerische und Angewandte Mathematik, BTU Cottbus-Senftenberg, Platz der Deutschen Einheit 1, D-03046 Cottbus, Germany, E-mail: Ludwig.Cromme@b-tu.de, Jens Kunath222Numerische und Angewandte Mathematik, BTU Cottbus-Senftenberg, Platz der Deutschen Einheit 1, D-03046 Cottbus, Germany and Andreas Krebs333Aerodynamik und Strömungslehre , BTU Cottbus-Senftenberg, Siemens-Halske-Ring 14, D-03046 Cottbus, Germany

Abstract

We present an algorithm to compute best least-squares approximations of discrete real-valued functions by first-degree splines (broken lines) with free knots. We demonstrate that the algorithm delivers after a finite number of steps a (global) best approximation. The analysis is complemented by remarks on programming and by a number of numerical examples including applications from medicine (MBC, MIC).

Key words: splines, first-degree splines, broken lines, splines with free knots, least-squares approximation, best approximation, minimal bactericidal concentration (MBC), minimal inhibitory concentration (MIC)

1 Introduction

Splines of a real variable with free knots form a highly nonlinear approximation family with varying dimension of the boundary (Cromme [2]). Even just local best approximations are therefore not easily computed (Cromme [3]) and globally convergent numerical methods are not known.

But for the special case of first-degree splines (broken lines) with free knots a globally convergent numerical method for the approximation of discrete data in the least-squares sense could be developed and is presented in the following sections.

The basic idea of the method is to reduce the complexity of the original problem: Instead of solving the original nonlinear problem on the whole data set a number of just linear problems on subsets (”segments”) is treated and assembled to a solution of the original problem.

The underlying idea is sketched in Fig. 1: In this figure the best $L_{2}$-approximation $s^{\ast}$ from the set of first-degree splines with $5$ free knots $t_{1},\ldots,t_{5}$ to discrete data $(x_{i},f_{i})$, $i=0,1,\ldots,18$, is plotted. Knots can be located between two data abscissae (such knots will be called interior knots in the following, in our example $t_{1}$ und $t_{4}$ are interior knots) or coincide with data abscissae (called data knots, in our example $t_{2}=x_{4},t_{3}=x_{10}$, and $t_{5}=x_{15}$ are data knots).

Basic is the observation that $s^{\ast}$ is not only best approximation on the whole data set, but also on subsets of the data called segments; for a formal introduction of notations see below. For example $\left.s^{\ast}\right|_{\left[x_{4},x_{12}\right]}$ is a best approximation from $S^{1}\left(x_{7},x_{10}\right)$, the first-degree splines with fixed knots $x_{7}$ and $x_{10}$, to the data $\left(x_{4},f_{4}\right),\ldots,\left(x_{12},f_{12}\right)$. For if a better approximation $\tilde{s}$ with fixed knots $x_{7}$ and $x_{10}$ to the data $\left(x_{4},f_{4}\right),\ldots,\left(x_{12},f_{12}\right)$ existed, $s^{\ast}$ could be modified in the considered segment in the direction of $\tilde{s}$. The interior knots $t_{1}$ and $t_{4}$ would be slightly moved in this process, the other knots would remain unchanged. The approximation error of the thus modified spline would be smaller than the approximation error of $s^{\ast}$ - contradicting the optimality of $s^{\ast}$ as best approximation.

So if these segments ranging from one interior knot to the next interior knot were known, we could determine $s^{\ast}$ by solving just linear (!) approximation problems with fixed knots segment for segment. But since the segments of a best approximation are not known in advance our algorithm examines systematically possible segmentations of the data abscissae to see if they might lead to a best approximation.

The implementation of this idea raises a number of questions which will be adressed after formal presentation of the algorithm in the following section.

Let us summarize the most important notations used: Let $\mathbb{P}_{m}$ denote the real-valued polynomials of degree smaller or equal $m$ and $[a,b]$ a real interval $a<b$. For $a=:t_{0}<t_{1}<\ldots,<t_{k}<t_{k+1}:=b$ and $m,k\in\mathbb{N}$ let

[TABLE]

be the set of splines of degree $m$ with $k$ fixed (simple) knots $t_{1},\ldots,t_{k}$. The splines of degree $m\in\mathbb{N}$ with at most $k\in\mathbb{N}$ free (simple) knots are the set

[TABLE]

$t_{i}$ is called an improper (or inactive) knot of a spline $s\in S^{m}_{k}[a,b]$ if $s$ is $m$-times differentiable in $t_{i}$, otherwise the knot is called proper (or active) knot.

Let $\mu+2$ real abscissae $a=:x_{0}<x_{1}<\ldots<x_{\mu+1}:=b$ and function values $f_{0},\ldots,f_{\mu+1}\in\mathbb{R}$ of a function $f:[a,b]\rightarrow\mathbb{R}$ be given. We define vectors $X:=\left(x_{0},\ldots,x_{\mu+1}\right)^{t}$, $F:=\left(f_{0},\ldots,f_{\mu+1}\right)^{t}$, and for $g:\left\{x_{0},\ldots,x_{\mu+1}\right\}\rightarrow\mathbb{R}$ the vector $g(X):=\left(g(x_{0}),\ldots,g(x_{\mu+1})\right)^{t}$. With the Euclidian norm $\left\|\cdot\right\|_{2}$ our problem can now be stated:

Definition 1.1

(Approximation problem) Find $s^{\ast}\in S^{1}_{k}[a,b]$ with

[TABLE]

In the following sections we generally assume $k\geq 1$ and $\mu\geq k+1\geq 2$, since otherwise the problem would be a purely polynomial approximation problem or the data could be reproduced exactly in $S_{k}^{1}[a,b]$.

2 The basic algorithm and convergence theorem

A spline $s\in S^{1}_{k}[a,b]$ partitions the data abscissae $x_{0}<x_{1}<\ldots<x_{\mu+1}$ naturally into segments each ranging from an interior knot (i.e. a knot situated between two data abscissae) to the next interior knot. Within a segment there can be more data knots, i.e. knots coinciding with data abscissae.

Such segmentations are basic for the numerical procedure presented below. In defining the segments it is advantageous to code the position of knots in a ”position vector”, i.e. the information wether a knot coincides with an abscissa and with which or between which abscissae a knot is located. To start we define position vectors for arbitrary first-degree splines. As will be seen later, the algorithm needs to deal only with ”regular” position vectors, see Defintion 2.6.

Definition 2.1

• (a)

A position vector $p$ of length $k$ to the data $x_{0},\ldots,x_{\mu+1}$ is a vector from $\mathbb{N}_{0}^{k}$ with

[TABLE]

[TABLE]

where $\mathbb{N}_{0}:=\mathbb{N}\cup\{0\}$ and [math] is considered even.

• (b)

For $s\in S^{1}_{k}[a,b]$, $a=:x_{0}<x_{1}<\ldots<x_{\mu+1}:=b$ and $x_{0}<t_{1}<t_{2}<\ldots<t_{k}<x_{\mu+1}$ the position vector $p\in\mathbb{N}_{0}^{k}$ of $s$ to the data $x_{0},\ldots,x_{\mu+1}$ is defined as

[TABLE]

Conclusion 2.2

• (a)

The position vector for $s^{\ast}$ from Fig. 1 is the vector $p=(6,13,19,24,29)$.

• (b)

A position vector of $s$ allows us to infer back the position of the knots of $s$: $t_{j}=x_{\frac{p(j)+1}{2}}$, if $p(j)$ is odd and $t_{j}\in\left(x_{\frac{p(j)}{2}},x_{\frac{p(j)}{2}+1}\right)$, if $p(j)$ is even. By Definition 2.1 (b) the position of knots is mapped one-to-one onto the set of position vectors.

• (c)

Successive components $p(j)$ and $p(j+1)$ of a position vector $p$ for $s\in S_{k}^{1}[a,b]$ are equal iff $t_{j}$ and $t_{j+1}$ lie in the same interval \big{(}x_{i},x_{i+1}\big{)}:

[TABLE]

• (d)

For every position vector $p$ according to Definition 2.1 (a) exists a spline $s\in S_{k}^{1}[a,b]$ with $p$ as position vector. To be more precice, there exists a whole class of splines all of which have $p$ as corresponding position vector.

Proof 1

We confine ourselves to the proof of (d). Let $p$ be a position vector according to Definition 2.1 (a). For odd $p(j)$ set $t_{j}:=x_{\frac{p(j)+1}{2}}$. For even $p(j)$ without less of generality let $p(j)$ to $p(j+q)$ be exactly the components of $p$ identical with $p(j)$, that is, $p(j)=p(j+1)=\ldots=p(j+q)$ and $p(j-1)<p(j)$, if $j\geq 2$ and $p(j+q)<p(j+q+1)$, if $j+q<k$. Choose $t_{j},\ldots,t_{j+q}$ such that

[TABLE]

holds. This implies $t_{1}<t_{2}<\ldots<t_{k}$ and each spline from $S_{k}^{1}[a,b]$ with these knots has position vector $p$. $\square$

The position of the knots (coded by the position vector) induces a segmentation of the data. Each segment ranges from one interior knot to the next, more precisely:

Definition/Conclusion 2.3

Let $p\in\mathbb{N}_{0}^{k}$ be a position vector for $s\in S^{1}_{k}[a,b]$ according to Definition 2.1:

Case 1:* All $p(j)$ are odd ($j=1,2,\ldots,k$). Then $s$ has no interior knots, all knots coalesce with date abscissae and we have*

[TABLE]

There is only one segment $S_{0}:=\left\{x_{0},x_{1},\ldots,x_{\mu+1}\right\}$.

• 2.

Case 2:* At least one component $p(j)$ of $p$ is even. Let $r\in\{1,\ldots,k\}$ be chosen maximal such that $p(j_{l})$ is even for $l=1,\ldots,r$ and suitable $j_{l}$ with $1\leq j_{1}<j_{2}<\ldots<j_{r}\leq k$. The knots $t_{j_{1}},\ldots,t_{j_{r}}$ are exactly the interior knots of $s$. Set*

[TABLE]

Let $r\geq 2$, $l\in\{1,\ldots,r-1\}$ hold. In case $i_{l}=i_{l+1}$, set $S_{l}:=\emptyset$. In case $i_{l}+1\leq i_{l+1}$ set

[TABLE]

If $i_{l}+1\leq i_{l+1}$ and in addition $q_{l}:=j_{l+1}-j_{l}-1>0$ holds then

[TABLE]

are the knots in segment $S_{l}$ coalescing with data abscissae. So the segments $S_{0}$ (to the left of first interior knot $t_{j_{1}}$) and $S_{r}$ (to the right of the last interior knot $t_{j_{r}}$) are left to be defined here (and also for $r=1$). $S_{0}$ consists of abscissae $x_{0},\ldots,x_{i_{1}}$ and contains - if $j_{1}-1>0$ - the knots

[TABLE]

Segment $S_{r}$ consists of abscissae $x_{i_{r}+1},\ldots,x_{\mu+1}$ and contains - if $j_{r}+1\leq k$ - the knots

[TABLE]

Example 2.4

Position vector $p=\left(6,1,3,19,24,29\right)$ for spline $s^{\ast}$ in Fig. 1 leads to segments $S_{0}=\left\{x_{0},x_{1},x_{2},x_{3}\right\}$ without knots, $S_{1}=\left\{x_{4},x_{5},\ldots,x_{12}\right\}$ with knots $t_{2}:=x_{7}$ and $t_{3}:=x_{10}$ and to segment $S_{2}=\left\{x_{13},x_{14},\ldots,x_{18}\right\}$ with the knot $t_{5}:=x_{15}$.

Remark 2.5

Because of Theorem 12 in [4] we can limit our search for a best approximation to splines with only simple knots and for which the additional features (a) to (d) from Theorem 12 in [4] hold; see also Appendix A. In Lemma 2.7 we will prove that these are exactly the splines with regular position vectors.

Definition 2.6

A position vector $p\in\mathbb{N}_{0}^{k}$ is called regular position vector of length $k$, if the following features hold:

• (a)

$\left\{p(1),\ldots,p(k)\right\}\subseteq\left\{1,2,\ldots,2\mu-1\right\}$*  and*

• (b)

\text{p(j) even}\quad\Rightarrow\quad\forall i\in\{1,\ldots,j-1\}:\;p(i)\neq p(j)-1,\;\forall i\in\{j+1,\ldots,k\}:\;p(i)\neq p(j)+1* .*

• (c)

For $p(j)$ and $p(j+l+1)$ ($l\geq 0$) both even and for $l\geq 1$ in addition $p(j+1),\ldots,p(j+l)$ odd the following inequality holds

[TABLE]

• (d)

\forall\,j\in\{1,\ldots,k-1\}:\;\;p(j+1)-p(j)\geq\left\{\begin{array}[]{lcl}4&,&\text{for p(j)$,$p(j+1) both even}\\ 2&,&\text{otherwise}\end{array}\right.**

Lemma 2.7

The position vector of a spline $s\in S^{k}[a,b]$ is regular iff $s$ exhibits features (a) thru (d) from Theorem 12 in [4].

Proof 2

Lemma 2.7 reflects the fact that properties (a) to (d) in Definition 2.6 for components $p(1),\ldots,p(k)$ guarantee features (a) to (d) in Theorem 12 in [4] and vice versa. As an example let us elaborate (d):

So let $p$ be a position vector of a spline $s\in S_{k}^{1}[a,b]$, $t_{j}$ and $t_{j+1}$ knots of $s$.

*Case 1: $t_{j}$ and *$t_{j+1}$are interior knots, i.e. $p(j)$ and $p(j+1)$ are even. Then at least two data abscissae are situated between $t_{j}$ and $t_{j+1}$ iff

[TABLE]

• -

Case 2: $t_{j}$ is an interior and $t_{j+1}$ a data knot, i.e., $p(j)$ is even and $p(j+1)$ is odd. Then at least two data abscissae lie between or on $t_{j}$ and $t_{j+1}$ iff

[TABLE]

The equivalence of feature (d) from Theorem 12 in **[4]** to the requirement

[TABLE]

is shown with similar arguments in

Case 3: $t_{j}$ is data knot, $t_{j+1}$ interior knot, that is, $p(j)$ is odd and $p(j+1)$ is even,

• -

Case 4: $t_{j}$ and $t_{j+1}$ are data knots, that is, $p(j)$ and $p(j+1)$ are odd.

Altogether this proves the equivalence of 2.6 (d) and Theorem 12 (d) from **[4]**. $\square$

Conclusion 2.8

For a regular position vector the following statements hold true:

[TABLE]

Proof 3

See Definition 2.6 (a) and (d).

Lemma 2.9

The number of regular position vectors grows exponentially in $\mu$.

Proof 4

According to Definition 2.6 a position vector is regular if all $k$ components $p(1)<p(2)<\ldots<p(k)$ are odd, i.e.

[TABLE]

In such cases all knots of related splines are data knots. So there are at least as many regular position vectors as there are subsets of $\left\{1,3,5,\ldots,2\mu-1\right\}$ with $k$ elements. According to standard results in combinatorics there are $\binom{\mu}{\mu-k}=\frac{\mu!}{(\mu-k)!k!}$ such subsets. With growing $\mu$ and $k=\frac{\mu}{2}$ for example exponential growth of the number of regular position vectors is therefore programmed. $\square$

The fast growth of the number of regular position vectors is evident in Table 1:

We can now specify the algorithm for the solution of Approximation problem 1.1: The algorithm examines systematically all regular position vectors to see if they might code the position of the knots of a best approximating spline $s\in S^{1}_{k}[a,b]$.

Algorithm 2.10

**

Set $\varrho:=\infty$.

• 2.

For the next regular position vector $p$ of length $k$ for the data abscissae $x_{0},\ldots,x_{\mu+1}$ which has not yet been examined carry out exactly one of the following two steps 3 and 4:

• 3.

*If all $p(i)$ for $i=1,2,\ldots,k$ are odd compute the best approximation *

$s\in S^{1}\left(x_{\frac{p(1)-1}{2}},\ldots,x_{\frac{p(k)-1}{2}}\right)$* to the data $\left(x_{0},f_{0}\right),\ldots,\left(x_{\mu+1},f_{\mu+1}\right)$. The spline $s$ is thus the solution of a discrete linear least-squares approximation problem with fixed knots. If \big{\|}F-s(X)\big{\|}_{2}<\varrho, set \varrho:=\big{\|}F-s(X)\big{\|}_{2}, $s^{\ast}:=s$. Continue with step 2.*

• 4.

If at least one component $p(i)$ is even, i.e., for a suitable $r\in\{1,\ldots,k\}$ there are exactly $r$ interior knots, carry out the following steps; notations are as in Definition 2.3, in addition $i_{0}:=-1,i_{r+1}:=\mu+2,j_{0}:=0,j_{r+1}:=k+1$, $q_{0}:=j_{1}-1$, $q_{r}:=k-j_{r}$:

• (a)

For each segment $S_{l}$ $\left(l=0,1,\ldots,r\right)$ compute the best least-squares approximation $s_{l}$ to the data $\left(x_{i_{l}+1},f_{i_{l}+1}\right),\ldots,\left(x_{i_{l+1}},f_{i_{l+1}}\right)$.   Here $s_{l}$ is to chosen from $\mathbb{P}_{1}$ if $q_{l}=0$ and

[TABLE]

if $q_{l}>0$.

• (b)

For all $l=1,\ldots,r$ check whether the splines $s_{l-1}$ and $s_{l}$ intersect in some point $z_{l}\in\left(x_{i_{l}},x_{i_{l}+1}\right)$. If this is the case for all $l$, $l=1,\ldots,r$ define

[TABLE]

If \big{\|}F-s(X)\big{\|}_{2}<\varrho set \varrho:=\big{\|}F-s(X)\big{\|}_{2}, $s^{\ast}:=s$ and procede with step 2. If for some $l\in\{1,\ldots,r\}$ splines $s_{l-1}$ and $s_{l}$ do not intersect in $\left(x_{i_{l}},x_{i_{l}+1}\right)$ reject the position vector $p$ under investigation and continue with step 2.

Procede with step 2.

• 5.

When all regular position vectors $p$ of length $k$ have been examined $s^{\ast}$ is the wanted best approximation and $\varrho$ the corresponding approximation error. End.

The proof that Algorithm 2.10 solves the approximation problem is based on Theorem 12 in [4]:

Theorem 2.11

Algorithm 2.10 stops after a finite number of steps and delivers with $s^{\ast}$ a best approximation

[TABLE]

Proof 5

Let $s^{\ast}$ be a best approximation whith the additional features from Theorem 12 in [4]. Because of Lemma 2.7 the position vector $p^{\ast}$ of $s^{\ast}$ is then regular and is therefore treated in Algorithm 2.10.

Furthermore, the following assertion holds:

Assertion 2.12

When Algorithm 2.10 elaborates on $p^{\ast}$ the function computed in step 3 or step 4 by the algorithm ist just $s^{\ast}$.

Proof 6

a) If step 3 is executed (all $p^{\ast}(i)$ are odd) then $s^{\ast}$ as best approximation from $S^{1}_{k}\left[a,b\right]$ is even more best approximation from the smaller set $S^{1}\left(x_{\frac{p(1)+1}{2}},\ldots,x_{\frac{p(k)+1}{2}}\right)$ which $s^{\ast}$ is also an element of.

b) If step 4 is executed every approximation $s_{i}$ on a segment is identical with $s^{\ast}$ in its domain. For simplicity we demonstrate this only for the first segment $S_{0}$.

Assumption: $s_{0}$ is a better linear approximation to the data $\left(x_{0},f_{0}\right),\ldots$, $\left(x_{\frac{p\left(j_{1}\right)}{2}},f_{\frac{p\left(j_{1}\right)}{2}}\right)$ than $s^{\ast}$. Any sufficiently small shift of $s^{\ast}$ in the direction of $s_{0}$ improves the approximation error. Because $\left(s^{\ast}\right)^{\prime}$ has a jump discontinuity in $t_{j_{1}}$ (see Theorem 12 (h) in [4]), the modification of $s^{\ast}$ in direction $s_{0}$ as described above results in a $s^{\ast}_{new}\in S_{k}^{1}[a,b]$ with $\left\|F-s^{\ast}_{neu}(X)\right\|_{2}<\left\|F-s^{\ast}(X)\right\|_{2}$ contradicting the optimality of $s^{\ast}$.

The above assumption must therefore be wrong and the assertion b) holds. $s^{\ast}$ is thus an element of $K\subset S^{1}_{k}[x_{0},x_{\mu+1}]$, the splines computed by the algorithm. Therefore the minimal approximation error is realized by the algorithm and a best approximation computed. $\square$

3 Realization, simplification and acceleration

As can be seen from Lemma 2.9 the computational cost and computing time grows exponentially with the number of data $\mu+2$ causing a corresponding memory demand. It ist advisable to restrain from computing and storing all (regular) position vectors since this would cause an exponential demand for memory. Instead, when a newly regular position vector is found, steps 3 und 4 should follow right away such that the position vector - once examined - can be deleted.

The algorithm is suitable for parallelization. Steps 3 and 4 can be executed for several regular position vectors in parallel.

Segments $S_{0}$ to $S_{r}$ and their best approximations $s_{0},\ldots,s_{r}$ could also be examined in parallel for intersection points. But this implementation is not advisable because it might lead to unnecessary computations: If the approxiomations $s_{0},s_{1},\ldots,s_{r}$ on segments $S_{0},S_{1},\ldots,S_{r}$ are computed in parallel and for example $s_{0}$ and $s_{1}$ do not intersect in \big{(}x_{i_{0}},x_{i_{0}+1}\big{)}, then the position vector under examination can already be dismissed and the computation of $s_{2},\ldots,s_{r}$ (for $r\geq 2$) in the parallel realization of the algorithm would be an additional and unnecessary effort.

In this respect more effective is sequential programming, starting step 4, computing $s_{0}$ and $s_{1}$, and then check whether they intersect in \big{(}x_{i_{0}},x_{i_{0}+1}\big{)}. In analogy a best approximation $s_{l+2}$ on the next segment $S_{l+2}$ is determined only if $s_{l}$ and $s_{l+1}$ have been found to intersect in \big{(}x_{i_{l}},x_{i_{l}+1}\big{)}. Otherwise, the position vector currently under investigation is sorted out and computing further approximations related to this position vector can be skipped.

Furthermore, the algorithm should be programmed such that if the program terminates prematurely because of time limitations at least the best approximation found so far is saved and handed over to the main program. If position vectors of good approximations are dealt with early the consequences of a premature termination can be further limited.

If there is prior knowledge about the position of knots of a best approximation related position vectors should be examined with priority.

Which clues for the optimal position of knots can be deduced from the data? Besides the divided differences of first order

[TABLE]

the (central) second order divided differences

[TABLE]

can be helpful here. The knots of a best-approximating first-order spline could be expected in the neighborhood of data abscissae $x_{j}$ where the $\Delta_{j}^{(1)}$ vary heavily and the absolute values of the $\Delta_{j}^{(2)}$ are ”large”. Related positions vectors should therefore be examined with priority.

4 Numerical examples and applications

In this section we present numerical examples to illustrate the performance of Algorithm 2.10. For each example we list the knots and the approximation error and present the solutions in graphical form, but refrain from enumerating the coefficients of the best approximations in all examples.

Example 4.1

Let $\mu=15$, $x_{i}:=i$ for $i=0,1,\ldots,16$, $f_{i}:=1$ for $i=0,1,\ldots,7,9,\ldots,16$ and for $f_{8}:=2$ be given. For $k=1,2,3,4,5$ Algorithm 2.10 computes best approximations $s_{k}\in S^{1}_{k}\left[t_{0},t_{16}\right]$ with knots $t_{1},\ldots,t_{k}$ given in the following table:

[TABLE]

For $k=1,2,3$ $s_{k}$ has exactly $k$ proper knots. To approximate the data for $k\geq 4$ only $3$ proper knots are needed and the remaining $k-3$ knots are not proper knots; see Fig. 3.

Example 4.2

*The so-called Titanium Heat Data have been dicused in literature and approximated from a number of approximation families; e.g. [5] and [6]. The data describe a property of titanium as function of temperature. The measuring points are equidistant in the interval $[595,1075]$. The least-squares approximation with broken lines delivers a good approximation as can seen in Fig. 4 and the following table:: *

[TABLE]

Example 4.3

We look at the equidistant data abscissae $x_{i}=-1+\frac{i}{10}$, $i=0,1,\ldots,20$ and the convex values $f_{i}:=x_{i}^{2}$, $i=0,1,\ldots,20$. Approximation of the data $(x_{i},f_{i})$ by a first-order spline with $k=5$ free knots preserves the convexity of the data as can be seen in Fig. 5.

Example 4.4

For $k=2,3$ the data $(x_{i},f_{i})$ from Fig. 6 were approximated by broken lines from $S^{1}_{k}\left[x_{0},x_{71}\right]$. The values $f_{i}$ are the aritmetic mean of measurements of forces measured in approximately equidistant time intervals when a sub-way switch ist set.

For $k=2$ a best approximation $s_{2}$ was computed with the (interior) knots $t_{1}=0.77277$ and $t_{2}=1.83478$ and approximation error $\left\|F-s_{2}(X)\right\|_{2}=686.36122$. Raising the number of free knots to $k=3$ lowers the approximation error to $\left\|F-s_{3}(X)\right\|_{2}=460.5584$ where the computed best approximation $s_{3}$ has the three (interior) knots $t_{1}=0.25851,t_{2}=0.87598$ and $t_{3}=1.81896$.

Example 4.5

Fig. 7 shows the best approximations from $S_{7}^{1}[x_{0},x_{19}]$ and from $S_{10}^{1}[x_{0},x_{19}]$ to $\mu+2=20$ data points $(x_{i},f_{i})$ with $x_{i}:=i+1$, $i=0,1,2,\ldots,20$. In computing the best approximation in $S_{7}^{1}[x_{0},x_{19}]$ Algorithm 2.10 inspected $6\,724\,520$ position vectors with $1\leq p(1)<p(2)<\ldots<p(10)\leq 2\mu-1$, only $795\,455$ of which are regular. For the approximation from $S_{10}^{1}[x_{0},x_{19}]$ the algorithm inspected $183\,579\,396$ position vectors with $1\leq p(1)<p(2)<\ldots<p(10)\leq 2\mu-1$, of which $1\,256\,465$ were regular.

Next, we present a medical application which originally motivated the development of Algorithm 2.10. The medical objective is to determine with an automated procedure the lowest concentration of an antibiotic necessary to kill (a defined high proportion of) the targeted bacteria, the so-called m**inimal bactericidal concentration (MBC). In clinical applications this helps to fight germs effectictely whith minimal unwanted side effects.

The measurements are determined as follows: An inital concentration $\kappa_{0}$ of the antibiotic is halved in each step. Germs are then exposed to each dilution $\kappa(j):=\kappa_{0}\cdot 2^{-j}$, $j=0,\ldots,z$, $z\in\mathbb{N}\setminus\{\,0\,\}$. The effect of concentration $\kappa(j)$ on the viability of the germs ist measured by a flourescence method which allows to derive the values

[TABLE]

for $j=0,\ldots,z$.

This gives us $z+1$ data points $\left(\kappa(0),\varphi_{0}\right),\ldots,\left(\kappa(z),\varphi_{z}\right)$. For better representation instead of concentration $\kappa(j)$ the logarithmic values $\log\left(\kappa(j)\right)$ are used, that ist, we make use of the relation $j\;=\;-\frac{\log(\kappa(j))-\log(k_{0})}{\log(2)}$ and plot the pairs

[TABLE]

See Fig. 8 for a typical graph of the data and a best approximating first-degree spline $s^{\ast}\in S^{1}_{2}\left[x_{0},x_{z}\right]$ with knots $t_{1}$ und $t_{2}$. The knot $t_{1}$ is interpreted from the medical perspective as minimal bactericidal concentration

[TABLE]

Also of medical interest is knot $t_{2}$ from which

[TABLE]

is derived, the m**inimal inhibitory concentration, the lowest concentration with noticeable bactericidal effect; see [1] for more details.

In the following examples the initial concentration was $\kappa_{0}:=256=2^{8}$, when $z=19$ and $\kappa_{0}:=128=2^{7}$, when $z=18$. The following numerical example shows typical data and a typical best approximation $s^{\ast}\in S^{1}_{2}\left[x_{0},x_{z}\right]$ with monotonically increasing $s^{\ast}$.

Example 4.6

With $x_{i}:=i$ for $i=0,1,\ldots,19$ the following medical values $f_{i}$ must be approximated:

\begin{array}[]{r|rrrrrrrrrrr}i&0&1&2&3&4&5&6&7&8&9&10\\ \hline\cr f_{i}&3.6273&3.381&3.0339&2.8414&2.7507&2.9006&2.941&2.9986&3.2127&3.8381&8.2629\end{array}**

\begin{array}[]{r|rrrrrrrrr}i&11&12&13&14&15&16&17&18&19\\ \hline\cr f_{i}&37.7363&84.0146&94.7914&98.7679&97.0424&98.0432&95.5602&99.0313&100\end{array}**

As best approximation from $S^{1}_{2}\left[x_{0},x_{19}\right]$ Algorithm 2.10 computed the broken line

[TABLE]

with the interior knots $t_{1}=10.2898$ and $t_{2}=12.2512$ and the approximation error $\left\|F-s^{\ast}(X)\right\|_{2}=5.7246$; see Fig. 8. This gives us the following values for MIC and MBC

[TABLE]

As seen in the above example the two knots $t_{1}$ and $t_{2}$ of a best approximation partition the medical data into three segments. In many cases the values are slightly monotonically increasing in the first and the last segment and ascending rapidly in the middle segment.

It may happen that no measurements are available from the middle segment. The following examples show that not always two (proper interior) knots are needed for a best approximation and the solutions need not to be unique.

Example 4.7

For $i=0,1,\ldots,18$ the medical data $\left(x_{i},f_{i}\right)$ are given by $x_{i}:=i$ and the values $f_{i}$ from the following table:

\begin{array}[]{r|rrrrrrrrrr}i&0&1&2&3&4&5&6&7&8&9\\ \hline\cr f_{i}&4.8245&5.0786&5.7781&6.105&5.9493&6.0516&5.589&5.5087&5.2563&4.5123\end{array}**

\begin{array}[]{r|rrrrrrrrr}i&10&11&12&13&14&15&16&17&18\\ \hline\cr f_{i}&97.8802&96.3044&95.6139&98.974&95.9425&96.0353&97.0482&98.5606&100\end{array}**

As best approximation from $S^{1}_{2}\left[x_{0},x_{18}\right]$ Algorithm 2.10 computed the first-degree spline

[TABLE]

with the data knots $t_{1}=9$ and $t_{2}=10$, see Fig. 9. The associated approximation error is $\left\|F-s^{\ast}(X)\right\|_{2}=4.24581$ and we compute $\text{MBC}=2^{7-t_{1}}=0.25$ and $\text{MIC}=2^{7-t_{2}}=0.125$.

Algorithm 2.10 computed just one best approximation, namely $s^{\ast}$, but there exist infinitely many best approximations to the given data: Each of the functions

[TABLE]

is a best approximation from $S^{1}_{2}\left[x_{0},x_{18}\right]$ to the data $\left(x_{0},f_{0}\right),\ldots,\left(x_{18},f_{18}\right)$ for arbitrary $y_{1},y_{2}$ with $t_{1}\leq y_{1}<y_{2}\leq t_{2}$. As limit case we get the discontinuous functions

[TABLE]

with jump discontinuity in $y$ for arbitrary $y\in\left(t_{1},t_{2}\right)$ as best approximations to the given data.

Example 4.8

For $i=0,1,\ldots,18$ we look at the data $(x_{i},f_{i})$ from Example 4.7 except that now $f_{9}:=7.5123$ and $f_{10}:=97.8802$. As best approximation from $S^{1}_{2}\left[x_{0},x_{18}\right]$ to the data pairs $\left(x_{0},f_{0}\right),\ldots,\left(x_{18},f_{18}\right)$ thus defined we get from Algorithm 2.10 the spline

[TABLE]

with the interior knot $t_{1}=8.98057$, the data knot $t_{2}=10$ and the minimal approximation error $\left\|F_{2}-s^{\ast}_{2}(X)\right\|_{2}=4.11872$, see Fig. 10. Here again, arbitrarily many solutions $s\left(y_{2},t\right)\in S^{1}_{2}\left[t_{0},t_{18}\right]$can be derived from $s^{\ast}$ for all of which

[TABLE]

must hold. This implies that for the knots $y_{1}$ and $y_{2}$ the inequality $t_{1}\leq y_{1}<x_{9}<y_{2}\leq t_{2}$ must be valid. Because of $s\left(y_{2},x_{9}\right)=f_{9}$ only one of the two knots can be freely chosen, the other one depends on that choice. If e.g. $y_{2}$ with $x_{9}<y_{2}<t_{2}$ ist chosen, then the splines

[TABLE]

also are best approximations, where $y_{1}$ ist determined as solution of $s_{1}(y_{2})=s_{2}(y_{2})$.

Measurement data not always meet the pattern of Example 4.6. Whether the interpretation of computed knots als MIC and MBC ist still a good idea must be decided from the medical perspective. In the next example we focus on such a situation.

Example 4.9

For $i=0,1,\ldots,19$ the medical data pairs $(x_{i},f_{i})$ are listed in the following table:

\begin{array}[]{r|rrrrrrrrrr}i&0&1&2&3&4&5&6&7&8&9\\ \hline\cr f_{i}&3.0354&3.1654&3.0862&3.0564&2.9804&2.9632&2.8198&3.1239&3.0576&2.9828\end{array}**

\begin{array}[]{r|rrrrrrrrrr}i&10&11&12&13&14&15&16&17&18&19\\ \hline\cr f_{i}&3.1498&3.5877&4.0296&6.6481&9.829&12.1237&30.1584&70.2245&89.7225&100\end{array}**

As best approximation from $S^{1}_{2}\left[t_{0},t_{19}\right]$ Algorithm 2.10 computed the spline

[TABLE]

witch the two (interior) knots $t_{1}=15.43646$ and $t_{2}=17.30953$ and the approximation error $\left\|F-s^{\ast}(X)\right\|_{2}=7.69589$, see Fig. 11. This leads to $\text{MBC}=2^{8-t_{1}}=0.00577$ and $\text{MIC}=2^{8-t_{2}}=0.00158$.

Appendix A: Theorem 12 from [4]

Theorem 12

Let $\mu\geq k+1\geq 2$ and $1\leq p\leq\infty$. For given data $(x_{0},f_{0}),\ldots$, $(x_{\mu+1},f_{\mu+1})$ with $x_{0}<x_{1}<\ldots<x_{\mu}<x_{\mu+1}$ and $f_{0},\ldots,f_{\mu+1}\in\mathbb{R}$ exists a best approximation $s^{\ast}\in S^{1}_{k}\left[x_{0},x_{\mu+1}\right]$

[TABLE]

with the following additional features where $t_{1},\ldots,t_{k}$ with $t_{0}:=x_{0}<t_{1}<\ldots<t_{k}<t_{k+1}:=x_{\mu+1}$ denote the knots of $s^{\ast}$:

• (a)

There are no knots in the boundary regions. More precisely:

[TABLE]

• (b)

Data abscissae neighboring to interior knots are not knots.

• (c)

Between two (not necessarily neighboring) interior knots of $s^{\ast}$ lie at least two data abscissae which are not knots of $s^{\ast}$.

• (d)

On or between neighboring knots lie at least two data abscissae:

[TABLE]

• (e)

If an interior knot is situated between two neighboring data abscissae, then no additional knot lies on or between these data abscissae. That is, the proposition ”$t_{j}$* is a knot of $s^{\ast}$ with $x_{q}<t_{j}<x_{q+1}$” implies $t_{j-1}<x_{q}$ and $t_{j+1}>x_{q+1}$.*

• (f)

Let $t_{j}$ be an interior knot. Then in each of the intervals $(-\infty,t_{j})$ and $(t_{j},\infty)$ there is a data abscissa from \big{\{}x_{1},\ldots,x_{\mu}\big{\}} which is not a knot.

• (g)

For $p<\infty$ we have: Between an interior and a neighboring data knot of $s^{\ast}$ lies either exactly one data abscissa $x_{q}$ which is then reproduced ($t_{j}<x_{q}<t_{j+1}$, $s^{\ast}(x_{q})=f_{q}$) or there exist at least two data abscissae $x_{q},x_{q+1}$ between the knots ($t_{j}<x_{q}<x_{q+1}<t_{j+1}$).

• (h)

All interior knots are proper knots, that is, the first derivative is discontinuous in all interior knots.