Riassumi il seguente testo. La struttura del riassunto sarà composta da un elenco di paragrafi di circa 100 parole ogniuno. Per ogni paragrafo aggiungi un titolo di massimo 4 o 5 parole che riassuma il contenuto del paragrafo.
Restituisci la risposta in formato JSON. La sintassi del JSON deve essere valida, quindi se nel titolo o nel testo devi inserire dei caratteri speciali, devi usare i caratteri escape, ad esempio \" oppure \'. La risposta deve includere soltanto il codice JSON per poterlo elaborare da un software. Non aggiungere nient'altro

Testo:
Nei giorni scorsi, Google ha annunciato un rebrand per la sua AI Bard che ora si chiama Google Gemini, riprendendo il nome del modello AI svelato lo scorso dicembre. Il rebrand riguarda tutte le versioni di Bard ed ora Gemini è disponibile anche con un’app per smartphone e tablet Android.
L’app non è ancora disponibile in Italia, ma arriverà certamente perché Google lo ha già annunciato. Non sarà più necessario passare dal browser per accedere all’assistente AI di Google. L’utilizzo dell’app di Gemini, però, ha una conseguenza: la disattivazione di Google Assistant. Gli utenti saranno, quindi, chiamati a una scelta.
Gemini disattiva Assistant
Come riportato da Android Authority, l’installazione dell’app di Gemini comporta una rinuncia a Google Assistant. Una volta installata l’applicazione dell’assistente AI, infatti, tutte le gesture utilizzate prima per Assistant saranno assegnate a Gemini. Gli utenti hanno la possibilità di ripristinare l’assistente vocale di Google, ma rinunciando a Gemini.
Da notare, però, riaprendo Gemini manualmente sarà necessario ripetere la procedura di configurazione. Parte di questa configurazione è proprio la sostituzione di Assistant. Di fatto, quindi, si crea un loop che obbliga l’utente disattivare l’assistente di Google ogni volta che intende utilizzare Gemini.
L’unica opzione alternativa per usare Gemini su Android, oltre all’app e alla conseguente disattivazione di Google Assistant, è rappresentata dal browser web, collegandosi all’indirizzo gemini.google.com. Il browser, almeno per il momento, è la soluzione migliore per usare l’AI della casa americana sul proprio smartphone Android senza dover rinunciare ad Assistant. L’accesso via browser è, inoltre, disponibile anche da altri dispositivi.
Come usare Gemini Advanced
Gemini Advanced è la versione più completa e potente di Gemini ed è accessibile tramite l’apposito piano in abbonamento, denominato Google One AI Premium. Questo piano ha un costo di 21,99 euro al mese (con i primi due mesi di prova gratuita) e offre anche 2 TB di spazio di archiviazione in cloud.
Come tutti i piani a pagamento di Google One, anche Google One AI Premium permette a chi lo sottoscrive di attivare un gruppo famiglia, con al massimo 5 membri.
La condivisione, però, non riguarda Gemini Advanced. La nuova tecnologia di Google, infatti, è disponibile solo per l’utente principale, quello che ha attiva l’abbonamento Google One AI Premium. Gli utenti invitati a condividere i vantaggi della sottoscrizione, invece, non potranno beneficiare della versione più potente di Gemini ma dovranno accontentarsi della versione standard.
Gli altri piani a pagamento di Google One non consentono l’accesso a Gemini Advanced, quindi se due o più persone all’interno dello stesso nucleo familiare vorranno sfruttare l’AI avanzata di Google, allora dovranno pagare tutti 21,99 euro al mese a testa.

JSON:
{
"paragrafi": [
{
"titolo": "Rebrand di Bard in Gemini",
"testo": "Google ha rinominato la sua AI Bard in Gemini, un potente assistente AI virtuale che può generare vari tipi di contenuti, disponibile ora anche tramite app per Android."
},
{
"titolo": "App Gemini e disattivazione di Assistant",
"testo": "L'app Gemini disattiva Google Assistant, questo costringe gli utenti a dover scegliere fra Gemini o Assistant."
},
{
"titolo": "Utilizzo di Gemini su Android",
"testo": "Gemini può essere utilizzato su Android tramite l'app (ancora non disponibile in italia, ma presto arriverà in base a quanto annunciato da Google) o dal browser web, se viene usata l'app Google Assistant verrà disabilitato."
},
{
"titolo": "Gemini Advanced",
"testo": "Gemini Advanced è la versione premium e più potente di Gemini, accessibile tramite l'abbonamento Google One AI Premium."
},
{
"titolo": "Condivisione di Gemini Advanced",
"testo": "Gemini Advanced è disponibile solo per l'utente principale dell'abbonamento Google One AI Premium e non può essere condiviso con i membri della famiglia. In tal caso sarà necessario pagare l'abbonamento per ogni membro della famiglia."
}
]
}

Testo:
Questo è il codice ufficiale di Stable Cascade. Forniamo script per l'addestramento e l'inferenza, oltre a una serie di modelli diversi che è possibile utilizzare.
Questo modello è costruito sull'architettura di Würstchen e la sua principale differenza rispetto ad altri modelli, come Stable Diffusion, è che lavora su uno spazio latente molto più piccolo. Perché è importante? Più piccolo è lo spazio latente, più veloce è l'inferenza e più economico è l'addestramento. Quanto è piccolo lo spazio latente? Stable Diffusion utilizza un fattore di compressione di 8, per cui un'immagine di 1024x1024 viene codificata a 128x128. Stable Cascade raggiunge un fattore di compressione di 42, il che significa che è possibile codificare un'immagine di 1024x1024 in 24x24, mantenendo ricostruzioni nitide. Il modello condizionale del testo viene quindi addestrato nello spazio latente altamente compresso. Le versioni precedenti di questa architettura hanno ottenuto una riduzione dei costi di 16 volte rispetto a Stable Diffusion 1.5.
Pertanto, questo tipo di modello è adatto per gli usi in cui l'efficienza è importante. Inoltre, con questo metodo sono possibili tutte le estensioni note come finetuning, LoRA, ControlNet, IP-Adapter, LCM ecc. Alcune di queste sono già state fornite (finetuning, ControlNet, LoRA) nelle sezioni dedicate all'addestramento e all'inferenza.
Inoltre, Stable Cascade ottiene risultati impressionanti, sia dal punto di vista visivo che da quello della valutazione. Secondo la nostra valutazione, Stable Cascade ottiene le migliori prestazioni sia nell'allineamento dei prompt che nella qualità estetica in quasi tutti i confronti. L'immagine qui sopra mostra i risultati di una valutazione umana che utilizza un mix di prompt parziali (link) ed estetici. In particolare, Stable Cascade (30 passi di inferenza) è stato confrontato con Playground v2 (50 passi di inferenza), SDXL (50 passi di inferenza), SDXL Turbo (1 passo di inferenza) e Würstchen v2 (30 passi di inferenza).

JSON:
{
"paragrafi": [
{
"titolo": "Stable Cascade",
"testo": "Stable Cascade è un modello di generazione di immagini che utilizza uno spazio latente più piccolo rispetto ad altri modelli, come Stable Diffusion, per un'inferenza più veloce e un addestramento più economico."
},
{
"titolo": "Spazio latente compresso",
"testo": "Stable Cascade utilizza un fattore di compressione di 42, che consente di codificare un'immagine di 1024x1024 in 24x24, mantenendo ricostruzioni nitide."
},
{
"titolo": "Efficienza e versatilità",
"testo": "Stable Cascade è adatto per gli usi in cui l'efficienza è importante e supporta estensioni come finetuning, LoRA, ControlNet e IP-Adapter."
},
{
"titolo": "Prestazioni impressionanti",
"testo": "Stable Cascade ottiene risultati eccellenti sia dal punto di vista visivo che della valutazione, superando altri modelli in termini di allineamento dei prompt e qualità estetica."
}
]
}

Testo:
ll between p(x) and f (x) except at the nodes.  Example 2.1. Consider these five data points: (0, 8), (1, 12), (3, 2), (4, 6), (8, 0). Construct and plot the interpolation polynomial using the two outermost points. Repeat this process by adding one additional point at a time until all the points are included. What conclusions can you draw? Solution: The first interpolation polynomial is the line between the outermost points (0, 8) and (8, 0). Then we added the points (3, 2), (4, 5), and (1, 12) in that order and plotted a curve for each additional point. All of these polynomials are shown in Figure 5. We were hoping for a smooth curve going through these points without wide fluctuations, but this did not happen. (Why?) It may seem counterintuitive, but as we added more points, the situation became worse instead of better! The reason for this comes from the nature of high-degree polynomials. A polynomial of degree n has n zeros. If all of these zero points are real, then the curve crosses the x-axis n times. The resulting curve must make many turns for this to happen, resulting in wild oscillations.  2.1. Dirichlet Functions. As a pathological example, consider the so-called Dirich-  let function f , defined to be 1 at each irrational point and 0 at each rational point. If we choose nodes that are rational numbers, then p(x) ≡ 0 and f (x) − p(x) = 0 for all rational values of x, but f (x) − p(x) = 1 for all irrational values of x.  However, if the function f is well-behaved, can we not assume that the differences |f (x) − p(x)| will be small when the number of interpolating nodes is large? The answer is still no, even for functions that possess continuous derivatives of all orders on the interval!  2.2. Runge Functions. A specific example of this remarkable phenomenon is pro-  vided by the Runge function:  (2.2.1)  f (x) = (1 + x2)−1  2. ERRORS IN POLYNOMIAL INTERPOLATION  Figure 6. Polynomial interpolant with nine equally spaced nodes.  on the interval [−5, 5]. Let pn be the polynomial that interpolates this function at n + 1 equally spaced points on the interval [−5, 5], including the endpoints. Then  lim n→∞  max −5≤x≤5  |f (x) − pn(x)| = ∞  Thus, the effect of requiring the agreement of f and pn at more and more points is to increase the error at nonnodal points, and the error actually increases beyond all bounds!  The moral of this example, then, is that polynomial interpolation of high degree with many nodes is a risky operation; the resulting polynomials may be very unsatisfactory as representations of functions unless the set of nodes is chosen with great care.  The reader can easily observe the phenomenon just described by using the pseu- docodes already developed in this chapter. In a more advanced study of this topic, it would be shown that the divergence of the polynomials can often be ascribed to the fact that the nodes are equally spaced. Again, contrary to intuition, equally distributed nodes are usually a very poor choice in interpolation. A much better choice for n + 1 nodes in [−1, 1] is the set of Chebyshev nodes: ã  ïÅ 2i + 1 2n + 2  ò  xi = cos  π  0 ≤ i ≤ n.  The corresponding set of nodes on an arbitrary interval [a, b] would be derived from a linear mapping to obtain  xi =  1 2  (a + b) +  1 2  (b − a) cos  ïÅ 2i + 1 2n + 2  ã  π  ò  0 ≤ i ≤ n.  Notice that these nodes are numbered from right to left. Since the theory does not depend on any particular ordering of the nodes, this is not troublesome.  A simple graph illustrates this phenomenon best. Again, consider Equation (2.2.1) on the interval [−5, 5]. First, we select nine equally spaced nodes and use routines Coef and Eval with an automatic plotter to graph p8. As shown in Figure 6, the resulting curve assumes negative values, which, of course, f (x) does not have! Adding more equally spaced nodes – and thereby obtaining a higher-degree polynomial?only makes matters worse with wilder oscillations. In Figure 7, nine Chebyshev nodes are used, and the resulting polynomial curve is smoother. However, cubic splines produce an even better curve fit.  The Chebyshev nodes are obtained by taking equally-spaced points on a semicircle  and projecting them down onto the horizontal axis, as in Figure 8.  157  158  CHAPTER 5.  INTERPOLATION  Figure 7. Polynomial interpolant with nine Chebyshev nodes.  Figure 8. Interpolation with Chebyshev points.  2.3. Theorems on Interpolation Errors. It is possible to assess the errors of interpolation by means of a formula that involves the (n + 1)st derivative of the function being interpolated. Here is the formal statement:  Theorem 2.2 (Interpolation errors I). If p is the polynomial of degree at most n that interpolates f at the n+1 distinct nodes x0, x1, . . . , xn belonging to an interval [a, b] and if f (n+1) is continuous, then for each x in [a, b], there is a ξ in (a, b) for which  (2.3.1)  f (x) − p(x) =  1 (n + 1)!  f (n+1)(ξ)  n (cid:89)  (x − xi)  i=0  Proof. Observe first that Equation (2.3.1) is obviously valid if x is one of the nodes xi because then both sides of the equation reduce to zero. If x is not a node, let it be fixed in the remainder of the discussion, and define  w(t) =  n (cid:89)  (t − xi)  (polynomial in the variable t)  (2.3.2)  i=0 f (x) − p(x) w(x) φ(t) = f (t) − p(t) − cw(t)  c =  (constant)  (function in the variable t)  2. ERRORS IN POLYNOMIAL INTERPOLATION  Observe that c is well defined because w(x) ̸= 0 (x is not a node). Note also that φ takes the value 0 at the n + 2 points x0, x1, . . ., xn, and x. Now invoke Rolle’s Theorem1, which states that between any two roots of φ, there must occur a root of φ′. Thus, φ′ has at least n + 1 roots. By similar reasoning, φ′′ has at least n roots, φ′′′ has at least n − 1 roots, and so on. Finally, it can be inferred that φ(n+1) must have at least one root. Let ξ be a root of φ(n+1). All the roots being counted in this argument are in (a, b). Thus,  0 = φ(n+1)(ξ) = f (n+1)(ξ) − p(n+1)(ξ) − cw(n+1)(ξ)  In this equation, p(n+1)(ξ) = 0 because p is a polynomial of degree ≤ n. Also, w(n+1)(ξ) = (n + 1)! because w(t) = tn+1 + (lower-order terms in t). Thus, we have  0 = f (n+1)(ξ) − c(n + 1)! = f (n+1)(ξ) −  (n + 1)! w(x)  [f (x) − p(x)]  This equation is a rearrangement of Equation (2.3.1).  A special case that often arises is the one in which the interpolation nodes are  equally spaced.  Lemma 2.3 (Upper Bound Lemma). Suppose that xi = a + ih for i = 0, 1, . . . , n  and that h = (b − a)/n. Then for any x ∈ [a, b]  (2.3.3)  n (cid:89)  i=0  |x − xi| ≤  1 4  hn+1n!  Proof. To establish this inequality, fix x and select j so that xj ≤ x ≤ xj+1. It is  an exercise in calculus to show that  (2.3.4)  |x − xj||x − xj+1| ≤  h2 4  Using Equation (2.3.4), we have  n (cid:89)  i=0  |x − xi| ≤  h2 4  j−1 (cid:89)  i=0  (x − xi)  n (cid:89)  i=j+2  (xi − x)  The sketch in Figure 9, showing a typical case of equally spaced nodes, may be helpful. Since xj ≤ x ≤ xj+1, we have further  n (cid:89)  i=0  |x − xi| ≤  h2 4  j−1 (cid:89)  i=0  (xj+1 − xi)  n (cid:89)  i=j+2  (xi − xj)  1Rolle’s Theorem: Let f be a function that is continuous on [a, b] and differentiable on (a, b). If  f (a) = f (b) = 0, then f ′(c) = 0 for some point c in (a, b).  159  □  160  CHAPTER 5.  INTERPOLATION  Figure 9. Typical location of x in equally spaced nodes.  Now use the fact that xi = a + ih. Then we have xj+1 − xi = (j − i + 1)h and xi − xj = (i − j)h. Therefore,  n (cid:89)  i=0  |x − xi| ≤  h2 4  hjhn−(j+2)+1  j−1 (cid:89)  i=0  (j − i + 1)  n (cid:89)  i=j+2  (i − j)  ≤  1 4  hn+1(j + 1)!(n − j)! ≤  1 4  hn+1n!  In the last step, we use the fact that if 0 ≤ j ≤ n − 1, then (j + 1)!(n − j)! ≤ n!. This, □ too, is left as an exercise. Hence, Inequality (2.3.3) is established.  We can now find a bound on the interpolation error. Theorem 2.4 (Interpolation errors II). Let f be a function such that f (n+1) is continuous on [a, b] and satisfies |f (n+1)(x)| ≤ M . Let p be the polynomial of degree ≤ n that interpolates f at n + 1 equally spaced nodes in [a, b], including the endpoints. Then on [a, b],  (2.3.5)  |f (x) − p(x)| ≤  1 4(n + 1)  M hn+1  where h = (b − a)/n is the spacing between nodes.  2.3.  Proof. Use Theorem 2.2 on interpolation errors and Inequality (2.3.3) in Lemma □  This theorem gives loose upper bounds on the interpolation error for different values of n. By other means, one can find tighter upper bounds for small values of n. If the nodes are not uniformly spaced then a better bound can be found by use of the Chebyshev nodes.  The error expression in polynomial interpolation can also be given in terms of divided  differences:  Theorem 2.5 (Interpolation errors III). If p is the polynomial of degree n that  interpolates the function f at nodes x0, x1, . . . , xn, then for any x that is not a node,  f (x) − p(x) = f [x0, x1, . . . , xn, x]  n (cid:89)  (x − xi)  i=0  Proof. Let t be any point, other than a node, where f (t) is defined. Let q be the polynomial of degree ≤ n + 1 that interpolates f at x0, x1, . . . , xn, t. By the Newton form of the interpolation formula [Equation (1.7.4)], we have  q(x) = p(x) + f [x0, x1, . . . , xn, t]  n (cid:89)  (x − xi)  i=0  2. PROBLEMS  Since q(t) = f (t), this yields at once  f (t) = p(t) + f [x0, x1, . . . , xn, t]  n (cid:89)  (t − xi)  i=0  The following theorem shows that there is a relationship between divided differences  and derivatives.  Theorem 2.6 (Divided differences and derivatives). If f (n) is continuous on [a, b]  and if x0, x1, . . . , xn are any n + 1 distinct points in [a, b], then for some ξ in (a, b),  f [x0, x1, . . . , xn] =  1 n!  f (n)(ξ)  Proof. Let p be the polynomial of degree ≤ n−1 that interpolates f at x0, x1, . . . , xn−1.  By Theorem 2.2 on interpolation errors, there is a point ξ such that  f (xn) − p(xn) =  1 n!  f (n)(ξ)  n−1 (cid:89)  i=0  (xn − xi)  By Theorem (2.5) on interpolation errors, we obtain  f (xn) − p(xn) = f [x0, x1, . . . , xn−1, xn]  n−1 (cid:89)  (xn − xi)  i=0  As an immediate consequence of this theorem, we observe that all high-order divided  differences are zero for a polynomial.  Corollary 2.7 (Divided Differences Corollary). If f is a polynomial of degree n,  then all of the divided differences f [x0, x1, . . . , xi] are zero for i ≥ n + 1.  Problems  (1) Use a divided-difference table to show that the following data can be repre-  sented by a polynomial of degree 3:  x −2 −1 4 1 y  3 0 11 16 13 −4  1  2  (2) For nonuniformly distributed nodes a = x0 < x1 < · · · < xn = b, where  h = max1≤i≤n (xi − xi−1), show that Inequality (2.3.3) is true.  (3) How accurately can we determine sin x by linear interpolation, given a table  of sin x to ten decimal places, for x in [0, 2] with h = 0.01?  (4) Let the function f (x) = ln x be approximated by an interpolation polynomial of degree 9 with ten nodes uniformly distributed in the interval [1, 2]. What bound can be placed on the error?  (5) Suppose cos x is to be approximated by an interpolating polynomial of degree n, using n + 1 equally spaced nodes in the interval [0, 1]. How accurate is the approximation? (Express your answer in terms of n.) How accurate is the approximation when n = 9? For what values of n is the error less than 10−7?  161  □  □  162  CHAPTER 5.  INTERPOLATION  (6) In interpolating with n + 1 equally spaced nodes on an interval, we could use xi = a + (2i + 1)h/2, where 0 ≤ i ≤ n − 1 and h = (b − a)/n. What bound can be given now for (cid:81)n i=0 |x − xi| when a ≤ x ≤ b? Note: We are not requiring the endpoints to be nodes.  (7) Does every polynomial p of degree at most n obey the following equation?  Explain why or why not.  p(x) =  n (cid:88)  p[x0, x1, . . . , xi]  i−1 (cid:89)  (x − xj)  i=0  j=0  Hint: Use the uniqueness of the interpolating polynomial.  Computer Problems  (1) Using 21 equally spaced nodes on the interval [−5, 5], find the interpolating polynomial p of degree 20 for the function f (x) = (x2 + 1)−1. Print the values of f (x) and p(x) at 41 equally spaced points, including the nodes. Observe the large discrepancy between f (x) and p(x).  (2) (Continuation) Perform the experiment in the preceding computer problem, using Chebyshev nodes xi = 5 cos(iπ/20), where 0 ≤ i ≤ 20, and nodes xi = 5 cos[(2i + 1)π/42], where 0 ≤ i ≤ 20. Record your conclusions.  (3) Let f (x) = max{0, 1 − x}. Sketch the function f . Then find interpolating polynomials p of degrees 2, 4, 8, 16, and 32 to f on the interval [−4, 4], using equally spaced nodes. Print out the discrepancy f (x) − p(x) at 128 equally spaced points. Then redo the problem using Chebyshev nodes.  (4) Why are the Chebyshev nodes generally better than equally spaced nodes in i=0(x−xi) that occurs  polynomial interpolation? The answer lies in the term (cid:81)n in the error formula. If xi = cos[(2i + 1)π/(2n + 2)], then  (cid:12) n (cid:12) (cid:89) (cid:12) (cid:12) (cid:12)  i=0  (x − xi)  (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)  ≤ 2−n  for all x ∈ [−1, 1]. Carry out a numerical experiment to test the given inequality for n = 3, 7, 15.  3. Approximations by Spline Functions  3.1. First-Degree and Second-Degree Splines. The history of spline functions is rooted in the work of draftsmen, who often needed to draw a gently turning curve between points on a drawing. This process is called fairing and can be accomplished with a number of ad hoc devices, such as the French curve, made of plastic and presenting a number of curves of different curvature for the draftsman to select. Long strips of wood were also used, being made to pass through the control points by weights laid on the draftsman’s table and attached to the strips. The weights were called ducks and the strips of wood were called splines, even as early as 1891. The elastic nature of the wooden strips allowed them to bend only a little while still passing through the prescribed points. The wood was, in effect, solving a differential equation and minimizing the strain energy. The latter is known to be a simple function of the curvature. The mathematical theory of these curves owes much to the early investigators, particularly Isaac Schoenberg in  3. APPROXIMATIONS BY SPLINE FUNCTIONS  Figure 10. First-degree spline function.  the 1940s and 1950s. Other important names associated with the early development of the subject (i.e., prior to 1964) are Garrett Birkhoff, C. de Boor, J. H. Ahlberg, E. N. Nilson, H. Garabedian, R. S. Johnson, F. Landis, A. Whitney, J. L. Walsh, and J. C. Holladay. The first book giving a systematic exposition of spline theory was the book by Ahlberg, Nilson, and Walsh [1967].  3.2. First-Degree Spline. A spline function is a function that consists of poly- nomial pieces joined together with certain smoothness conditions. A simple example is the polygonal function (or spline of degree 1), whose pieces are linear polynomials joined together to achieve continuity, as in Figure 9.2. The points t0, t1, . . . , tn at which the function changes its character are termed knots in the theory of splines. Thus, the spline function shown in Figure 10 has eight knots. Such a function appears somewhat complicated when defined in explicit terms. We are forced to write  (3.2.1)  S(x) =       S0(x) S1(x) ...  x ∈ [t0, t1] x ∈ [t1, t2]  ... Sn−1(x) x ∈ [tn−1, tn]  where  (3.2.2)  Si(x) = aix + bi  because each piece of S(x) is a linear polynomial. Such a function S(x) is piecewise lin- ear. If the knots t0, t1, . . . , tn were given and if the coefficients a0, b0, a1, b1, . . . , an−1, bn−1 were all known, then the evaluation of S(x) at a specific x would proceed by first deter- mining the interval that contains x and then using the appropriate linear function for that interval.  If the function S defined by Equation (3.2.1) is continuous, we call it a first-degree  spline. It is characterized by the following three properties.  Definition 3.1 (Spline of Degree 1). A function S is called a spline of degree 1 if: (1) The domain of S is an interval [a, b]. (2) S is continuous on [a, b]. (3) There is a partitioning of the interval a = t0 < t1 < · · · < tn = b such that S  is a linear polynomial on each subinterval [ti, ti+1].  163  164  CHAPTER 5.  INTERPOLATION  Outside the interval [a, b], S(x) is usually defined to be the same function on the left of a as it is on the leftmost subinterval [t0, t1] and the same on the right of b as it is on the rightmost subinterval [tn−1, tn], namely, S(x) = S0(x) when x < a and S(x) = Sn−1(x) when x > b.  Continuity of a function f at a point s can be defined by the condition  f (x) = lim x→s− Here, limx→s+ means that the limit is taken over x values that converge to s from above s; that is, (x − s) is positive for all x values. Similarly, limx→s− means that the x values converge to s from below.  lim x→s+  f (x) = f (s)  Example 3.2. Determine whether this function is a first-degree spline function:  S(x) =       x ∈ [−1, 0] x 1 − x x ∈ (0, 1) 2x − 2 x ∈ [1, 2]  Solution: The function is obviously piecewise linear but is not a spline of degree 1 because it is discontinuous at x = 0. Notice that limx→0+ S(x) = limx→0(1 − x) = 1, whereas limx→0− S(x) = limx→0 x = 0.  The spline functions of degree 1 can be used for interpolation. Suppose the following  table of function values is given:  t1 x t0 y y0 y1  · · · · ·  tn tn  There is no loss of generality in supposing that t0 < t1 < . . . < tn because this is only a matter of labeling the knots.  The table can be represented by a set of n+1 points in the plane, (t0, y0), (t1, y1),. . ., (tn, yn), and these points have distinct abscissas. Therefore, we can draw a polygonal line through the points without ever drawing a vertical segment. This polygonal line is the graph of a function, and this function is obviously a spline of degree 1. What are the equations of the individual line segments that make up this graph?  By referring to Figure 11 and using the point-slope form of a line, we obtain  (3.2.3)  Si(x) = yi + mi(x − ti)  on the interval [ti, ti+1], where mi is the slope of the line and is therefore given by the formula  mi =  yi+1 − yi ti+1 − ti  Notice that the function S that we are creating has 2n parameters in it: the n coef- ficients ai and the n constants bi in Equation (3.2.2). On the other hand, exactly 2n conditions are being imposed, since each constituent function Si must interpolate the data at the ends of its subinterval. Thus, the number of parameters equals the number of conditions. For the higher-degree splines, we shall encounter a mismatch in these two numbers; the spline of degree k will have k − 1 free parameters for us to use as we wish in the problem of interpolating at the knots.  3. APPROXIMATIONS BY SPLINE FUNCTIONS  Figure 11. First-degree spline: linear Si(x).  The form of Equation (3.2.3) is better than that of Equation (3.2.2) for the practical evaluation of S(x) because some of the quantitiesx − ti must be computed in any case simply to determine which subinterval contains x. If t0 ≤ x ≤ tn then the interval [ti, ti+1] containing x is characterized by the fact that x − ti is the first of the quantities x − tn−1, x − tn−2, . . . , x − t0 that is nonnegative.  Procedure 35 is a function procedure that utilizes n + 1 table values (ti, yi) in linear arrays (ti) and (yi), assuming that a = t0 < t1 < . . . < tn = b. Given an x value, the routine returns S(x) using Equations (3.2.1) and (3.2.3). If x < t0, then S(x) = y0 + m0(x − t0); if x > tn, then S(x) = yn−1 + mn−1(x − tn−1) .  Algorithm 35  function Spline1(n, (ti), (yi), x) for i = n − 1 : −1 : 0 do if x − ti ≥ 0 then  break  end if  end for Spline1 ← yi + (x − ti)[(yi+1 − yi)/(ti+1 − ti)]  end function  3.3. Modulus of Continuity. To assess the goodness of fit when we interpolate a function with a first-degree spline, it is useful to have something called the modulus of continuity of a function f . Suppose f is defined on an interval [a, b]. The modulus of continuity of f is  ω(f ; h) = sup{|f (u) − f (v)| : a ≤ u ≤ v ≤ b, |u − v| ≤ h}  Here, sup is the supremum, which is the least upper bound of the given set of real numbers. The quantity ω(f ; h) measures how much f can change over a small interval  165  166  CHAPTER 5.  INTERPOLATION  of width h. If f is continuous on [a, b], then it is uniformly continuous, and ω(f ; h) will tend to zero as h tends to zero. If f is not continuous, ω(f ; h) will not tend to zero. If f is differentiable on (a, b) (in addition to being continuous on [a, b]) and if f ′(x) is bounded on (a, b), then the Mean Value Theorem can be used to get an estimate of the modulus of continuity: If u and v are as described in the definition of ω(f ; h), then  |f (u) − f (v)| = |f ′(c)(u − v)| ≤ M1|u − v| ≤ M1h Here, M1 denotes the maximum of |f ′(x)| as x runs over (a, b). For example, if f (x) = x3 and [a, b] = [1, 4], then we find that ω(f ; h) ≤ 48h.  Theorem 3.3 (First-degree Polynomial Accuracy Theorem). If p is the first-degree polynomial that interpolates a function f at the endpoints of an interval [a, b], then with h = b − a, we have  |f (x) − p(x)| ≤ ω(f ; h) a ≤ x ≤ b  Proof. The linear function p is given explicitly by the formula Å x − a b − a  ã  Å b − x b − a  ã  p(x) =  f (b) +  f (a)  Hence,  f (x) − p(x) =  Å x − a b − a  ã  [f (x) − f (b)] +  Å b − x b − a  ã  [f (x) − f (a)]  Then we have  Å x − a b − a Å x − a b − a ïÅ x − a b − a = ω(f ; h)  |f (x) − p(x)| ≤  ≤  =  ã  Å b − x b − a ã  ã  |f (x) − f (b)| +  |f (x) − f (a)|  ã  Å b − x b − a ãò  ω(f ; h) +  ω(f ; h)  ã  Å b − x b − a  +  ω(f ; h)  From this basic result, one can easily prove the following one, simply by applying  the basic inequality to each subinterval.  Theorem 3.4 (First-degree Spline Accuracy Theorem). Let p be a first-degree spline having knots a = x0 < x1 < · · · < xn = b. If p interpolates a function f at these knots, then with h = maxi(xi − xi−1), we have  |f (x) − p(x)| ≤ ω(f ; h) a ≤ x ≤ b  If f ′ or f ′′ exist and are continuous, then more can be said, namely,  |f (x) − p(x)| ≤ M1  |f (x) − p(x)| ≤ M2  h 2 h2 8  a ≤ x ≤ b  a ≤ x ≤ b  In these estimates, M1 is the maximum value of |f ′(x)| on the interval, and M2 is the maximum of |f ′′(x)|.  □  3. APPROXIMATIONS BY SPLINE FUNCTIONS  The first theorem tells us that if more knots are inserted in such a way that the maximum spacing h goes to zero, then the corresponding first-degree spline will converge uniformly to f . Recall that this type of result is conspicuously lacking in the polynomial interpolation theory. In that situation, raising the degree and making the nodes fill up the interval will not necessarily ensure that convergence takes place for an arbitrary continuous function.  3.4. Second-Degree Splines. Splines of degree higher than 1 are more compli- cated. We now take up the quadratic splines. Let’s use the letter Q to remind ourselves that we are considering piecewise quadratic functions. A function Q is a second-degree spline if it has the following properties.  Definition 3.5 (Spline of Degree 2). A function Q is called a spline of degree 2  if:  (1) The domain of Q is an interval [a, b]. (2) Q and Q′ are continuous on [a, b]. (3) There are points ti (called knots) such that a = t0 < t1 < . . . < tn = b and Q  is a polynomial of degree at most 2 on each subinterval [ti, ti+1].  In brief, a quadratic spline is a continuously differentiable piecewise quadratic func- tion, where quadratic includes all linear combinations of the basic functions x (cid:55)→ 1, x, x2.  Example 3.6. Determine whether the following function is a quadratic spline:  Q(x) =       x2 −10 ≤ x ≤ 0 −x2 0 ≤ x ≤ 1 1 − 2x 1 ≤ x ≤ 20  Solution: The function is obviously piecewise quadratic. Whether Q and Q′ are contin- uous at the interior knots can be determined as follows:  lim x→0− lim x→1− lim x→0− lim x→1−  Q(x) = lim x→0− Q(x) = lim x→1− Q′(x) = lim x→0− Q′(x) = lim x→1−  x2 = 0  (−x2) = −1  2x = 0  (−2x) = −2  lim x→0+ lim x→1+ lim x→0+ lim x→1+  Q(x) = lim x→0+ Q(x) = lim x→1+ Q′(x) = lim x→0+ Q′(x) = lim x→1+  (−x2) = 0  (1 − 2x) = −1  (−2x) = 0  (−2) = −2  Consequently, Q(x) is a quadratic spline.  3.5. Interpolating Quadratic Spline Q(x). Quadratic splines are not used in applications as often as are natural cubic splines, which are developed in the next sec- tion. However, the derivations of interpolating quadratic and cubic splines are similar enough that an understanding of the simpler second-degree spline theory will allow one to grasp easily the more complicated third-degree spline theory. We want to emphasize that quadratic splines are rarely used for interpolation, and the discussion here is pro- vided only as preparation for the study of higher-order splines, which are used in many applications.  167  168  CHAPTER 5.  INTERPOLATION  Proceeding now to the interpolation problem, suppose that a table of values has  been given:  x t0 t2 y y0 y1 y2  t1  · · · · ·  tn tn  We shall assume that the points t0, t1, . . . , tn, which we think of as the nodes for the interpolation problem, are also the knots for the spline function to be constructed. Later, another quadratic spline interpolant is discussed in which the nodes for interpolation are different from the knots.  A quadratic spline, as just described, consists of n separate quadratic functions x (cid:55)→ aix2 + bix + ci, one for each subinterval created by the n + 1 knots. Thus, we start with 3n coefficients. On each subinterval [ti, ti+1], the quadratic spline function Qi must satisfy the interpolation conditions Qi(ti) = yi and Qi(ti+1) = yi+1. Since there are n such subintervals, this imposes 2n conditions. The continuity of Q does not add any additional conditions. (Why?) However, the continuity of Q′ at each of the interior knots gives n − 1 more conditions. Thus, we have 2n + n − 1 = 3n − 1 conditions, or one condition short of the 3n conditions required. There are a variety of ways to impose this additional condition; for example, Q′(t0) = 0 or Q′′(t0) = 0.  We now derive the equations for the interpolating quadratic spline, Q(x). The value of Q′(t0) is prescribed as the additional condition. We seek a piecewise quadratic function  (3.5.1)  Q(x) =       Q0(x) Q1(x) ... Qn−1(x)  t0 ≤ x ≤ t1 t1 ≤ x ≤ t2  tn−1 ≤ x ≤ tn  which is continuously differentiable on the entire interval [t0, tn] and which interpolates the table; that is, Q(ti) = yi for 0 ≤ i ≤ n.  Since Q′ is continuous, we can put zi := Q′(ti). At present, we do not know the  correct values of zi; nevertheless, the following must be the formula for Qi:  (3.5.2)  Qi(x) =  zi+1 − zi 2(ti+1 − ti)  (x − ti)2 + zi(x − ti) + yi  To see that this is correct, verify that Qi(ti) = yi, Q′ i(ti+1) = zi+1. These three conditions define the function Qi uniquely on [ti, ti+1] as given in Equation (3.5.2).  i(ti) = zi, and Q′  Now, for the quadratic spline function Q to be continuous and to interpolate the table of data, it is necessary and sufficient that Qi(ti+1) = yi+1 for i = 0, 1, . . . , n − 1 in Equation (3.5.2). When this equation is written out in detail and simplified, the result is  Å yi+1 − yi ti+1 − ti This equation can be used to obtain the vector [z0, z1, . . . , zn]⊤, starting with an arbi- trary value for z0. We summarize with an algorithm: Algorithm: Quadratic Spline Interpolation at the Knots  ã  (3.5.3)  zi+1 = −zi + 2  0 ≤ i ≤ n − 1  3. APPROXIMATIONS BY SPLINE FUNCTIONS  Figure 12. Subbotin quadratic splines (t0 = τ0, t3 = τ4).  (1) Determine [z0, z1, . . . , zn]⊤ by selecting z0 arbitrarily and computing z1, z2, . . . , zn  recursively by Formula (3.5.3).  (2) The quadratic spline interpolating function Q is given by Formulas (3.5.1) and  (3.5.2).  3.6. Subbotin Quadratic Spline. A useful approximation process, first proposed by Subbotin [1967], consists of interpolation with quadratic splines, where the nodes for interpolation are chosen to be the first and last knots and the midpoints between the knots. Remember that knots are defined as the points where the spline function is permitted to change in form from one polynomial to another. The nodes are the points where values of the spline are specified. In the Subbotin quadratic spline function, there are n + 2 interpolation conditions and 2(n − 1) conditions from the continuity of Q and Q′. Hence, we have the exact number of conditions needed, 3n, to define the quadratic spline function completely.  We outline the theory here, leaving details for the reader to fill in. Suppose that  knots a = t0 < t1 < . . . < tn = b have been specified; let the nodes be the points  ® τ0 = t0 τi = 1  2 (ti + ti−1)  τn+1 = tn 1 ≤ i ≤ n  We seek a quadratic spline function Q that has the given knots and takes prescribed values at the nodes:  Q(τi) = yi,  0 ≤ i ≤ n + 1  as in Figure 12. The knots create n subintervals, and in each of them, Q can be a different quadratic polynomial. Let us say that on [ti, ti+1], Q is equal to the quadratic polynomial Qi. Since Q is a quadratic spline, it and its first derivative should be continuous. Thus, zi ≡ Q′(ti) is well defined, although as yet we do not know its values. It is easy to see that on [ti, ti+1], our quadratic polynomial can be represented in the form  (3.6.1)  Qi(x) = yi+1 +  1 2  (zi+1 + zi)(x − τi+1) +  1 2hi  (zi+1 − zi)(x − τi+1)2  in which hi = ti+1 − ti. To verify the correctness of Equation (3.6.1), we must check that Qi(τi+1) = yi+1, Q′ i(ti+1) = zi+1. When the polynomial pieces  i(ti) = zi, and Q′  169  170  CHAPTER 5.  INTERPOLATION  Q0, Q1, . . . , Qn−1 are joined together to form Q, the result may be discontinuous. Hence, we impose continuity conditions at the interior knots:  lim x→t− i  Qi−1(x) = lim x→t+ i  Qi(x)  1 ≤ i ≤ n − 1.  The reader should carry out this analysis, which leads to  (3.6.2)  hi−1zi−1 + 3(hi−1 + hi)zi + hizi+1 = 8(yi+1 − yi)  1 ≤ i ≤ n − 1  The first and last interpolation conditions must also be imposed:  Q(τ0) = y0 Q(τn+1) = yn+1  These two equations lead to  3h0z0 + h0z1 = 8(y1 − y0)  hn−1zn−1 + 3hn−1zn = 8(yn+1 − yn) The system of equations governing the vector z = [z0, z1, . . . , zn]⊤ then can be written in the matrix form              3h0 h0  h0 3(h0 + h1) h1  h1 3(h1 + h2) . . .  h2 . . . hn−2  . . . 3(hn−2 + hn−1)  hn−1 hn−1  3hn−1                          z0 z1 z2 ... zn−1 zn              = 8              y1 − y0 y2 − y1 y3 − y2 ... yn − yn−1 yn+1 − yn  This system of n + 1 equations in n + 1 unknowns can be conveniently solved by procedure T ri. After the z vector has been obtained, values of Q(x) can be computed from Equation (3.6.1).  Problems  (1) Determine whether this function is a first-degree spline:  S(x) =       x 0.5 + 2(x − 0.5) x + 1.5  −1 ≤ x ≤ 0.5 0.5 ≤ x ≤ 2 2 ≤ x ≤ 4  (2) Let t0 < t1 < . . . < tn. Construct first-degree spline functions G0, G1, . . . , Gn by requiring that Gi vanish at t0, t1, . . . , ti−1, ti+1, . . . , tn but that Gi(ti) = 1. Show that the first-degree spline function that interpolates f at t0, t1, . . . , tn is (cid:80)n  i=0 f (ti)Gi(x).  (3) Prove that the derivative of a quadratic spline is a first-degree spline. (4) Define f (x) = 0 if x < 0 and f (x) = x2 if x ≥ 0. Show that f and f ′ are continuous. Show that any quadratic spline with knots t0, t1, . . . , tn is of the form  ax2 + bx + c +  n−1 (cid:88)  dif (x − ti)  i=1 (5) What equations must be solved if a quadratic spline function Q that has knots 2 (ti + ti+1) for  t0, t1, . . . , tn is required to take prescribed values at points 1 0 ≤ i ≤ n − 1?              3. COMPUTER PROBLEMS  (6) Are these functions quadratic splines? Explain why or why not.  ® 0.1x2  0 ≤ x ≤ 1 1 ≤ x ≤ 1.3  (a) Q(x) =  9.3x2 − 18.4x + 9.2 ® −x2 −100 ≤ x ≤ 0  (b) Q(x) =  x  0 ≤ x ≤ 100     x −50 ≤ x ≤ 1 x2 1 ≤ x ≤ 2 2 ≤ x ≤ 50 4 (7) Is S(x) = |x| a first-degree spline? Why or why not?  (c) Q(x) =    Computer Problems  (1) Rewrite procedure Spline1 so that a binary search is used to find the desired interval. Test the revised code. What are the advantages and/or disadvantages of a binary search compared to the procedure in the notes? A binary search is similar to the bisection method in that we choose tk with k = (i + j)/2 or k = (i + j + 1)/2 and determine whether is in [ti, tk] or [tk, tj].  (2) A piecewise bilinear polynomial that interpolates points (x, y) specified in  a rectangular grid is given by  p(x, y) =  (lijzi+1,j+1 + li+1,j+1zij) − (li+1,jzi,j+1 + li,j+1zi+1,j) (xi+1 − xi)(yj+1 − yj)  where lij = (xi − x)(yj − y). Here xi ≤ x ≤ xi+1 and yj ≤ y ≤ yj+1. The given grid (xi, yj) is specified by strictly increasing arrays (xi) and (yj) of length n and m, respectively. The given values zij at the grid points (xi, yj) are contained in the n × m array (zij), shown in Figure 13. Write  real function Bi Linear((xi), n, (yj), m, (zij), x, y)  to compute the value of p(x, y). Test this routine on a set of 5 × 10 unequally spaced data points. Evaluate Bi Linear at four grid points and five nongrid points.  Figure 13  171  172  CHAPTER 5.  INTERPOLATION  (3) Write an adaptive spline interpolation procedure. The input should be a func- tion f , an interval [a, b], and a tolerance ε. The output should be a set of knots a = t0 < t1 < . . . < tn = b and a set of function values yi = f (ti) such that the first-degree spline interpolating function S satisfies |S(x) − f (x)| ≤ ε whenever x is any point xij = ti + j(ti+1 − tj)/10 for 0 ≤ i ≤ n − 1 and 0 ≤ j ≤ 9.  4. Natural Cubic Splines  4.1. Introduction. The first- and second-degree splines discussed in the preceding section, though useful in certain applications, suffer an obvious imperfection: Their low- order derivatives are discontinuous. In the case of the first-degree spline (or polygonal line), this lack of smoothness is immediately evident because the slope of the spline may change abruptly from one value to another at each knot. For the quadratic spline, the discontinuity is in the second derivative and is therefore not so evident. But the curvature of the quadratic spline changes abruptly at each knot, and the curve may not be pleasing to the eye.  The general definition of spline functions of arbitrary degree is as follows. Definition 4.1 (Spline of Degree k). A function S is called a spline of degree k  if:  (1) The domain of S is an interval [a, b]. (2) S, S′, S′′, . . . , S(k−1) are all continuous functions on [a, b]. (3) There are points ti (the knots of S) such that a = t0 < t1 < . . . < tn = b and  such that S is a polynomial of degree at most k on each subinterval [ti, ti+1].  Observe that no mention has been made of interpolation in the definition of a spline function. Indeed, splines are such versatile functions that they have many applications other than interpolation.  Higher-degree splines are used whenever more smoothness is needed in the approx- imating function. From the definition of a spline function of degree k, we see that such a function will be continuous and have continuous derivatives S′, S′′, . . ., S(k−1). If we want the approximating spline to have a continuous mth derivative, a spline of degree at least m+1 is selected. To see why, consider a situation in which knots t0 < t1 < . . . < tn have been prescribed. Suppose that a piecewise polynomial of degree m is to be defined, with its pieces joined at the knots in such a way that the resulting spline S has m contin- uous derivatives. At a typical interior knot t, we have the following circumstances: To the left of t, S(x) = p(x); to the right of t, S(x) = q(x), where p and q are mth-degree polynomials. The continuity of the mth derivative S(m) implies the continuity of the lower-order derivatives S(m−1), S(m−2), . . ., S′, S. Therefore, at the knot t,  lim x→t− from which we conclude that  S(k)(x) = lim x→t+  S(k)(x)  0 ≤ k ≤ m  (4.1.1)  lim x→t−  p(k)(x) = lim x→t+  q(k)(x)  0 ≤ k ≤ m  Since p and q are polynomials, their derivatives of all orders are continuous, and so Equation (4.1.1) is the same as  p(k)(t) = q(k)(t)  0 ≤ k ≤ m  4. NATURAL CUBIC SPLINES  This condition forces p and q to be the same polynomial because by Taylor’s Theorem,  p(x) =  m (cid:88)  k=0  1 k!  p(k)(t)(x − t)k =  m (cid:88)  k=0  1 k!  q(k)(t)(x − t)k = q(x)  This argument can be applied at each of the interior knots t1, t2, . . . , tn−1, and we see that S is simply one polynomial throughout the entire interval from t0 to tn. Thus, we need a piecewise polynomial of degree m + 1 with at most m continuous derivatives to have a spline function that is not just a single polynomial throughout the entire interval. (We already know that ordinary polynomials usually do not serve well in curve fitting.) The choice of degree most frequently made for a spline function is 3. The resulting splines are termed cubic splines. In this case, we join cubic polynomials together in such a way that the resulting spline function has two continuous derivatives everywhere. At each knot, three continuity conditions will be imposed. Since S, S′, and S′′ are continuous, the graph of the function will appear smooth to the eye. Discontinuities, of course, will occur in the third derivative but cannot be easily detected visually, which is one reason for choosing degree 3. Experience has shown, moreover, that using splines of degree greater than 3 seldom yields any advantage. For technical reasons, odd-degree splines behave better than even-degree splines (when interpolating at the knots). Finally, a very elegant theorem, to be proved later, shows that in a certain precise sense, the cubic interpolating spline function is the best interpolating function available. Thus, our emphasis on the cubic splines is well justified.  4.2. Natural Cubic Spline. We turn next to interpolating a given table of func- tion values by a cubic spline whose knots coincide with the values of the independent variable in the table. As earlier, we start with the table:  t1 x t0 y y0 y1  tn . . . . . . yn  The ti’s are the knots and are assumed to be arranged in ascending order.  The function S that we wish to construct consists of n cubic polynomial pieces:  S(x) =       S0(x) S1(x) ... Sn−1(x)  t0 ≤ x ≤ t1 t1 ≤ x ≤ t2  ...  tn−1 ≤ x ≤ tn  In this formula, Si denotes the cubic polynomial that will be used on the subinterval [ti, ti+1]. The interpolation conditions are  S(ti) = yi  0 ≤ i ≤ n  The continuity conditions are imposed only at the interior knots t1, t2, . . . , tn−1. (Why?) These conditions are written as  lim x→t− i  S(k)(ti) = lim x→t+ i  S(k)(ti) k = 0, 1, 2  173  174  CHAPTER 5.  INTERPOLATION  It turns out that two more conditions must be imposed to use all the degrees of freedom available. The choice that we make for these two extra conditions is  S′′(t0) = S′′(tn) = 0 The resulting spline function is then termed a natural cubic spline. Additional ways to close the system of equations for the spline coefficients are periodic cubic splines and clamped cubic splines. A clamped spline is a spline curve whose slope is fixed at both end points: S′(t0) = d0 and S′(tn) = dn. A periodic cubic spline has S(t0) = S(tn), S′(t0) = S′(tn), and S′′(t0) = S′′(tn). For all continuous differential functions, clamped and natural cubic splines yield the least oscillations about the function f that it interpolates.  (4.2.1)  We now verify that the number of conditions imposed equals the number of coef- ficients available. There are n + 1 knots and hence n subintervals. On each of these subintervals, we shall have a different cubic polynomial. Since a cubic polynomial has four coefficients, a total of 4n coefficients are available. As for the conditions imposed, we have specified that within each interval the interpolating polynomial must go through two points, which gives 2n conditions. The continuity adds no additional conditions. The first and second derivatives must be continuous at the n − 1 interior points, for 2(n − 1) more conditions. The second derivatives must vanish at the two endpoints for a total of 2n + 2(n − 1) + 2 = 4n conditions.  Example 4.2. Derive the equations of the natural cubic interpolating spline for the  following table:  x −1 0 1 y  1 2 −1  Solution: Our approach is to determine the parameters a, b, c, d, e, f, g, and h so that S(x) is a natural cubic spline, where  S(x) =  ® S0(s) = ax3 + bx2 + cx + d  x ∈ [−1, 0]  S1(s) = ex3 + f x2 + gx + h x ∈ [0, 1]  where the two cubic polynomials are S0(x) and S1(x). From these interpolation con- ditions, we have interpolation conditions S(−1) = S0(−1) = −a + b − c + d = 1, S(0) = S0(0) = d = 2, S(0) = S1(0) = h = 2, and S(1) = S1(1) = e + f + g + h = −1. Taking the first derivatives, we obtain ® S′ 0(x) = 3ax2 + 2bx + c 1(x) = 3ex2 + 2f x + g S′ 0(0) = S′  S′(x) =  From the continuity condition of S′, we have S′ taking the second derivatives, we obtain ® S′′ S′′  1(0), and we set c = g. Next  0 (x) = 6ax + 2b 1 (x) = 6ex + 2f  S′′(x) =  From the continuity condition of S′′, we have S′′ 1 (0), and we let b = f . For S to be a natural cubic spline, we must have S′′ 1 (1) = 0, and we obtain 3a = b and 3e = −f . From all of these equations, we obtain a = −1, b = −3, c = −1, d = 2, e = 1, f = −3, g = −1, and h = 2.  0 (0) = S′′ 0 (−1) = 0 and S′′  4. NATURAL CUBIC SPLINES  4.3. Algorithm for Natural Cubic Spline. From the previous example, it is evident that we need to develop a systematic procedure for determining the formula for a natural cubic spline, given a table of interpolation values. This is our objective in the material on the next several pages.  Since S′′ is continuous, the numbers  zi := S′′(ti)  0 ≤ i ≤ n  are unambiguously defined. We do not yet know the values z1, z2, . . . , zn−1, but, of course, z0 = zn = 0 by Equation (4.2.1).  If the zi’s were known, we could construct S as now described. On the interval [ti, ti+1], S′′ is a linear polynomial that takes the values zi and zi+1 at the endpoints. Thus,  (4.3.1)  S′′  i (x) =  zi+1 hi  (x − ti) +  zi hi  (ti+1 − x)  with hi = ti+1 − ti for 0 ≤ i ≤ n − 1. To verify that Equation (4.3.1) is correct, notice that S′′ is linear in x. If this is integrated twice, we obtain Si itself:  i (ti) = zi, S′′  i (ti+1) = zi+1, and S′′ i  Si(x) =  zi+1 6hi  (x − ti)3 +  zi 6hi  (ti+1 − x)3 + cx + d  where c and d are constants of integration. By adjusting the integration constants, we obtain a form for Si that is easier to work with, namely,  (4.3.2)  Si(x) =  zi+1 6hi  (x − ti)3 +  zi 6hi  (ti+1 − x)3 + Ci(x − ti) + Di(ti+1 − x)  where Ci and Di are constants. If we differentiate Equation (4.3.2) twice, we obtain Equation (4.3.1).  The interpolation conditions Si(ti) = yi and Si(ti+1) = yi+1 can be imposed now to determine the appropriate values of Ci and Di. The reader should do so and verify that the result is  (4.3.3)  Si(x) =  zi+1 6hi Å yi+1 hi  zi 6hi ã  (x − ti)3 +  hi 6  +  −  zi+1  (ti+1 − x)3  (x − ti) +  Å yi hi  −  hi 6  zi  ã  (ti+1 − x)  When the values z0, z1, . . . , zn have been determined, the spline function S(x) is ob- tained from equations of this form for S0(x), S1(x), . . . , Sn−1(x). We now show how to determine the zi’s. One condition remains to be imposed – namely, the continuity of S′. At the interior knots ti for 1 ≤ i ≤ n − 1, we must have S′ i(ti), as can be seen in Figure 14.  i−1(ti) = S′  We have, from Equation (4.3.3),  S′  i(x) =  zi+1 2hi  (x − ti)2 −  zi 2hi  (ti+1 − x)2 +  yi+1 hi  −  hi 6  zi+1 −  yi hi  +  hi 6  zi  This gives  (4.3.4)  S′  i(ti) = −  hi 6  zi+1 −  hi 3  zi + bi  175  176  CHAPTER 5.  INTERPOLATION  Figure 14. Cubic spline: adjacent pieces Si−1 and Si.  where  (4.3.5)  bi =  1 hi  (yi+1 − yi)  Analogously, we have  S′  i−1(ti) =  hi−1 6  zi−1 +  hi−1 3  zi + bi−1  When these are set equal to each other, the resulting equation can be rearranged as  hi−1zi−1 + 2(hi−1 + hi)zi + hizi+1 = 6(bi − bi−1)  for 1 ≤ i ≤ n − 1. By letting  (4.3.6)  ui = 2(hi−1 + hi) vi = 6(bi − bi−1)  we obtain a tridiagonal system of equations:  (4.3.7)       z0 = 0 hi−1zi−1 + uizi + hizi+1 = vi zn = 0  1 ≤ i ≤ n − 1  to be solved for the zi’s. The simplicity of the first and last equations is a result of the natural cubic spline conditions S′′(t0) = S′′(tn) = 0. Now consider System (4.3.7) in matrix form:             1 0 h0 u1 h1 h1 u2 . . .  h2 . . . . . . hn−2 un−1 hn−1 0  1                        z0 z1 z2 ... zn−1 zn             =             0 v1 v2 ... vn−1 0             4. NATURAL CUBIC SPLINES  On eliminating the first and last equations, we have     u1 h1 h1 u2 . . .  h2 . . . . . . hn−3 un−2 hn−2 hn−2 un−1         (4.3.8)                  z1 z2 ... zn−2 zn−1           =           v1 v2 ... vn−2 vn−1           which is a symmetric tridiagonal system of order n − 1. We could use procedure T ri to solve this system. However, we can design an algorithm specifically for it. In Gaussian elimination without pivoting, the forward elimination phase would modify the ui’s and vi’s as follows:  (cid:40)  h2 i−1 ui−1  ui ← ui − vi ← vi − hi−1vi−1  ui−1  i = 2, 3, . . . , n − 1  The back substitution phase yields  ® zn−1 ← vn−1 un−1 zi ← vi−hizi+1  ui  i = n − 2, n − 3, . . . , 1  Putting all this together leads to the following algorithm, designed especially for the tridiagonal System (4.3.8). Algorithm: Solving the Natural Cubic Spline Tridiagonal System Directly Given the interpolation points (ti, yi) for i = 0, 1, . . . , n:  (1) Compute for i = 0, 1, . . . , n − 1:  ® hi = ti+1 − ti  bi = 1 hi  (yi+1 − yi)  (2) Set  ® u1 = 2(h0 + h1) v1 = 6(b1 − b0)  and compute inductively for i = 2, 3, . . . , n − 1:  (cid:40)  h2 i−1 ui−1  ui = 2(hi + hi−1) − vi = 6(bi − bi−1) − hi−1vi−1  ui−1  (3) Set  ® zn = 0 z0 = 0  and compute inductively for i = n − 1, n − 2, . . . , 1:  zi =  vi − hizi+1 ui  This algorithm conceivably could fail because of divisions by zero in steps 2 and 3.  177  178  CHAPTER 5.  INTERPOLATION  Therefore, let us prove that ui  ̸= 0 for all i.  It is clear that u1 > h1 > 0.  ui−1 > hi−1, then ui > hi because  ui = 2(hi + hi−1) −  h2 i−1 ui−1  > 2(hi + hi−1) − hi−1 > hi  Then by induction, ui > 0 for i = 1, 2, . . . , n − 1.  Equation (4.3.3) is not the best computational form for evaluating the cubic poly-  nomial Si(x). We would prefer to have it in the form  (4.3.9)  Si(x) = Ai + Bi(x − ti) + Ci(x − ti)2 + Di(x − ti)3  because nested multiplication can then be utilized.  Notice that Equation (4.3.9) is the Taylor expansion of Si about the point ti. Hence,  Ai = Si(ti), Bi = S′  i(ti), Ci =  1 2  S′′  i (ti), Di =  1 6  S′′′  i (ti)  Therefore, Ai = yi and Ci = zi/2. The coefficient of x3 in Equation (4.3.9) is Di, whereas the coefficient of x3 in Equation (4.3.3) is (zi+1 − zi)/6hi. Therefore,  Di =  1 6hi  (zi+1 − zi)  Finally, Equation (4.3.4) provides the value of S′  i(ti), which is  Bi = −  hi 6  zi+1 −  hi 3  zi +  1 hi  (yi+1 − yi)  Thus, the nested form of Si(x) is  (4.3.10)  Si(x) = yi + (x − ti)  Å  Bi + (x − ti)  Å zi 2  +  1 6hi  (x − ti)(zi+1 − zi)  ãã  4.4. Pseudocode for Natural Cubic Spline. We now write routines for deter- mining a natural cubic spline based on a table of values and for evaluating this function at a given value. First, we use Algorithm 36 for directly solving the tridiagonal System (4.3.8). This procedure, called Spline3 Coef , takes n + 1 table values (ti, yi) in arrays (ti) and (yi) and computes the zi’s, storing them in array (zi). Intermediate (working) arrays (hi), (bi), (ui), and (vi) are needed.  Now a procedure called Spline3 Eval is written for evaluating Equation (4.3.10), the natural cubic spline function S(x), for x a given value. The procedure Spline3 Eval first determines the interval [ti, ti+1] that contains x and then evaluates Si(x) using the nested form of this cubic polynomial: The function Spline3 Eval can be used repeatedly with different values of x after one call to procedure Spline3 Coef . For example, this would be the procedure when plotting a natural cubic spline curve. Since procedure Spline3 Coef stores the solution of the tridiagonal system corresponding to a particular spline function in the array (zi), the arguments n, (ti), (yi), and (zi) must not be altered between repeated uses of Spline3 Eval.  If  4. NATURAL CUBIC SPLINES  Algorithm 36  function Spline3 Coef(n, (ti), (yi), (zi))  for i = 0 : n − 1 do hi ← ti+1 − ti bi ← (yi+1 − yi)/hi  end for u1 ← 2(h0 + h1) v1 ← 6(b1 − b0) for i = 2 : n − 1 do  ui ← 2(hi + hi−1) − h2 vi ← 6(bi − bi−1) − hi−1vi−1/ui−1  i−1/ui−1  end for zn ← 0 for i = n − 1 : −1 : 1 do zi ← (vi − hizi+1)/ui  end for z0 ← 0 end function  Algorithm 37  function Spline3 Eval(n, (ti), (yi), (zi), x)  for i = n − 1 : −1 : 0 do if x − ti ≥ 0 then  break  end if  end for h ← ti+1 − ti tmp ← (zi/2) + (x − ti)(zi+1 − zi)/(6h) tmp ← −(h/6)(zi+1 + 2zi) + (yi+1 − yi)/h + (x − ti)(tmp) Spline3 Eval ← yi + (x − ti)(tmp)  end function  4.5. Smoothness Property. Why do spline functions serve the needs of data fit- ting better than ordinary polynomials? To answer this, one should understand that interpolation by polynomials of high degree is often unsatisfactory because polynomials may exhibit wild oscillations. Polynomials are smooth in the technical sense of possess- ing continuous derivatives of all orders, whereas in this sense, spline functions are not smooth.  Wild oscillations in a function can be attributed to its derivatives being very large. Consider the function whose graph is shown in Figure 15. The slope of the chord that joins the points p and q is very large in magnitude. By the Mean-Value Theorem, the slope of that chord is the value of the derivative at some point between p and q. Thus, the derivative must attain large values. Indeed, somewhere on the curve between p and q, there is a point where f ′(x) is large and negative. Similarly, between q and r, there is a point where f ′(x) is large and positive. Hence, there is a point on the curve between  179  180  CHAPTER 5.  INTERPOLATION  Figure 15. Wildly oscillating function.  p and r where f ′′(x) is large. This reasoning can be continued to higher derivatives if there are more oscillations. This is the behavior that spline functions do not exhibit. In fact, the following result shows that from a certain point of view, natural cubic splines are the best functions to use for curve fitting.  Theorem 4.3 (Cubic Spline Smoothness Theorem). If S is the natural cubic spline function that interpolates a twice-continuously differentiable function f at knots a = t0 < t1 < . . . < tn = b, then  (cid:90) b  [S′′(x)]2 dx ≤  (cid:90) b  [f ′′(x)]2 dx  a Proof. To verify the assertion about [S′′(x)]2, we let  a  g(x) = f (x) − S(x)  so that g(ti) = 0 for 0 ≤ i ≤ n, and  f ′′ = S′′ + g′′  Now  (cid:90) b  (f ′′)2 dx =  (cid:90) b  (S′′)2 dx +  (cid:90) b  (g′′)2 dx + 2  (cid:90) b  S′′g′′ dx  a If the last integral were 0, we would be finished because then  a  a  a  (cid:90) b  (f ′′)2 dx =  (cid:90) b  (S′′)2 dx +  (cid:90) b  (g′′)2 dx′ ≥  (cid:90) b  (S′′)2 dx  a  a  a  a  We apply the technique of integration by parts to the integral in question to show that it is 0. We have  (cid:90) b  a  S′′g′′ dx = S′′g′(cid:12) (cid:12) (cid:12)  b  a  −  (cid:90) b  a  S′′′g′ dx = −  (cid:90) b  a  S′′′g′ dx  4. PROBLEMS  Here, use has been made of the fact that S is a natural cubic spline; that is, S′′(a) = 0 and S′′(b) = 0. Continuing, we have  (cid:90) b  S′′′g′ dx =  n−1 (cid:88)  (cid:90) ti+1  S′′′g′ dx  a  i=0  ti  Since S is a cubic polynomial in each interval [ti, ti+1], its third derivative there is a constant, say ci. So  (cid:90) b  S′′′g′ dx =  n−1 (cid:88)  ci  (cid:90) ti+1  g′ dx =  n−1 (cid:88)  ci[g(ti+1) − g(ti)] = 0  a  i=0  ti  i=0  because g vanishes at every knot.  The interpretation of the integral inequality in the theorem is that the average value of [S′′(x)]2 on the interval [a, b] is never larger than the average value of this expression with any twice-continuous function f that agrees with S at the knots. The quantity [f ′′(x)]2 is closely related to the curvature of the function f .  Problems  (1) Determine the parameters a, b, c, d, and e such that S is a natural cubic spline: ® a + b(x − 1) + c(x − 1)2 + d(x − 1)3 x ∈ [0, 1] x ∈ [1, 2]  S(x) =  (x − 1)3 + ex2 − 1  (2) Determine the values of a, b, c, and d such that f is a cubic spline and such  that (cid:82) 2  0 [f ′′(x)]2 dx is a minimum: ® 3 + x − 9x3  f (x) =  a + b(x − 1) + c(x − 1)2 + d(x − 1)3  0 ≤ x ≤ 1 1 ≤ x ≤ 2  (3) Suppose S(x) is an mth-degree interpolating spline function over the interval  [a, b] with n + 1 knots a = t0 < t1 < . . . < tn = b. (a) How many conditions are needed to define S(x) uniquely over [a, b]? (b) How many conditions are defined by the interpolation conditions at the  knots?  (c) How many conditions are defined by the continuity of the derivatives? (d) How many additional conditions are needed so that the total equals the  number in part (a)?  (4) Determine the coefficients in the function  S(x) =  ® x3 − 1  −9 ≤ x ≤ 0  ax3 + bx2 + cx + d 0 ≤ x ≤ 5  such that it is a cubic spline that takes the value 2 when x = 1.  (5) Assume that a = x0 < x1 < . . . < xm = b. Describe the function f that interpolates a table of values (xi, yi), where 0 ≤ i ≤ m, and that minimizes the expression (cid:82) b  a |f ′(x)| dx.  181  □  182  CHAPTER 5.  INTERPOLATION  (6) Let knots t0 < t1 < . . . < tn, and let numbers yi and zi be given. Determine formulas for a piecewise cubic function f that has the given knots such that f ′′(x) = f (ti) = yi (0 ≤ i ≤ n), limx→t+ zi (1 ≤ i ≤ n). Why is f not generally a cubic spline?  f ′′(x) = zi (0 ≤ i ≤ n−1), and limx→t−  i  i  (7) A periodic cubic spline having knots t0, t1, . . . , tn is defined as a cubic spline function S(x) such that S(t0) = S(tn), S′(t0) = S′(tn), and S′′(t0) = S′′(tn). It would be used to fit data that are known to be periodic. Carry out the analysis necessary to obtain a periodic cubic spline interpolant for the table  x t0 t1 y y0 y1  · · · · ·  tn yn  assuming that yn = y0.  (8) Given a differentiable function f and knots t0 < t1 < · · · < tn, show how to obtain a cubic spline S that interpolates f at the knots and satisfies the end conditions S′(t0) = f ′(t0) and S′(tn) = f ′(tn). Note: This procedure produces a better fit to f when applicable. If f ′ is not known, finite-difference approximations to f ′(t0) and f ′(tn) can be used.  Computer Problems  (1) Rewrite and test procedure Spline3 Coef using procedure T ri. Use the sym-  metry of the (n − 1) × (n − 1) tridiagonal system.  (2) Let S be the cubic spline function that interpolates f (x) = (x2 + 1)−1 at 41 equally spaced knots in the interval [−5, 5]. Evaluate S(x)−f (x) at 101 equally spaced points on the interval [0, 5].  (3) Draw a free-form curve on graph paper, making certain that the curve is the graph of a function. Then read values of your function at a reasonable number of points, say, 10 − 50, and compute the cubic spline function that takes those values. Compare the freely drawn curve to the graph of the cubic spline.  (4) Write a program to estimate (cid:82) b  a f (x) dx, assuming that we know the values of f at only certain prescribed knots a = t0 < t1 < · · · < tn = b. Approximate f first by an interpolating cubic spline, and then compute the integral of it using Equation (4.3.3).  (5) Write a procedure to estimate f ′(x) for any x in [a, b], assuming that we know  only the values of f at knots a = t0 < t1 < · · · < tn = b.  CHAPTER 6  Introduction to Numerical Optimization  1. Linear Least Squares Problems[2, Chap.6]  This chapter is concerned with methods for solving the algebraic problem  min x  ∥b − Ax∥2,  where the dimensions of the real matrix and vectors are  A ∈ Rm×n, x ∈ Rn, b ∈ Rm, m ≥ n.  We assume that A has full column rank. Notice that in the overdetermined case, m > n, there is typically no x satisfying Ax = b exactly, even in the absence of roundoff error. The least squares problem arises often in many diverse application fields, especially where data fitting is required. Instances arise in machine learning, computer vision, and computer graphics applications, to name but a few. In computer vision people may want to match local invariant features of cluttered images under arbitrary rotations, scalings, change of brightness and contrast, and so on. Such actions are parametrized and the best values for the parameters are then sought to match the data. In computer graphics a parametrization of a surface mesh in three dimensions may be sought that yields a morphing of one animal into another, say. In such applications the question of just how to parametrize, or how to generate an efficient predictive model A in the present terminology, is often the crux of the matter and is far from trivial. But once the method of parametrization is determined, what follows is an instance of the problem considered here. Other applications seek to find an approximating function v(t, x) depending on a continuous variable t that fits data pairs (ti, bi), i = 1, . . . , m.  2. Least Squares and the Normal Equations  The space spanned by the columns of our m × n matrix A (i.e., all vectors z of the form z = Ay) is generally of dimension at most n, and we further assume that the dimension equals n, so that A has full column rank. In other words, we assume that its columns are linearly independent. Notice also that in the overdetermined case, m > n, b generally does not lie in the range space of A.  2.1. Deriving the normal equations. Let us rewrite the problem that we aim  to solve as  min x  1 2  ∥b − Ax∥2.  We have squared the normed expression, thus getting rid of the square root sign, and multiplied it by 1/2: this will not change the minimizer, in the sense that the same solution coefficients xj will be obtained. Notice that we have dropped the subscript 2  183  184  CHAPTER 6.  INTRODUCTION TO NUMERICAL OPTIMIZATION  Figure 1. Matrices and vectors and their dimensions in l2 data fitting.  in the norm notation: there is no other norm here to get confused by. Finally, we define the residual vector as usual by  r = b − Ax.  Writing these matrix and vectors explicitly we have  x =       x1 ... xn      , A =       a11 ... am1  · ·  · ·  a1n ... amn      , b =       b1 ... bm      , r =       r1 ... rm      .  Note that the matrix A is m × n, with m > n and perhaps m ≫ n, so it is “long and skinny”, and correspondingly we do not expect r to vanish at the optimum. See Figure 1.  We have a minimization problem for a smooth scalar function in several variables,  given by  min x  ψ(x), where ψ(x) =  1 2  ∥r∥2.  2. LEAST SQUARES AND THE NORMAL EQUATIONS  The necessary conditions for a minimum are obtained by setting the derivatives of ψ with respect to each unknown xk to zero, yielding  ∂ ∂xk  ψ(x) = 0,  k = 1, . . . , n.  Since  ψ(x) =  1 2  ∥r∥2 =  1 2  m (cid:88)  i=1  Ñ  bi −  n (cid:88)  j=1  aijxj  é2  ,  the conditions for a minimum yield  ∂ ∂xk  ψ(x) =  m (cid:88)  i=1      Ñ  bi −  n (cid:88)  j=1  aijxj  é  (−aik)     = 0  for k = 1, 2, . . . , n. The latter expression can be rewritten as  m (cid:88)  aik  n (cid:88)  aijxj =  m (cid:88)  aikbi  k = 1, . . . , n.  i=1  j=1  i=1  In matrix-vector form this expression looks much simpler; it reads  A⊤Ax = A⊤b.  This system of n linear equations in n unknowns is called the normal equations. Note that B = A⊤A can be much smaller in size than A; see Figure 1. The matrix B is symmetric positive definite given that A has full column rank.  2.2. Least Squares Solution Uniqueness. Is the solution of the normal equa- tions really a minimizer (and not, say, a maximizer) of the least squares norm? and if yes, is it the global minimum? The Least Squares Theorem given on the next page says it is. The answer is affirmative because the matrix B is positive definite.  To see this (completing the proof of the theorem), note first that our objective  function is a quadratic in n variables. Indeed, we can write  ψ(x) =  1 2  (b − Ax)⊤(b − Ax) =  1 2  x⊤Bx − b⊤Ax +  1 2  ∥b∥2.  Theorem 2.1 (Least Squares). The least squares problem  min x  ∥Ax − b∥2,  where A has full column rank, has a unique solution that satisfies the normal equations  Ä A⊤A  ä  x = A⊤b.  185  186  CHAPTER 6.  INTRODUCTION TO NUMERICAL OPTIMIZATION  Figure 2. Discrete least squares approximation.  2.3. Solving via the Normal Equations. It is important to realize that x, which has been our vector argument above, is now specified as the solution for the normal equations, which is indeed the solution of the least squares minimization problem.  Furthermore, if A has full column rank, then B = A⊤A is symmetric positive definite. The least squares problem has therefore been reduced, at least in principle, to that of solving an n × n system of linear equations.  The beauty does not end here, geometrically speaking. Notice that for the corre-  sponding residual at optimum we have  A⊤(b − Ax) = A⊤r = 0.  Hence we seek a solution satisfying that the residual is orthogonal to the column space of A. Since a picture is better than a thousand words (or at least a thousand bytes), Figure 2 is provided to illustrate the projection.  Our solution is given by x = (A⊤A)−1A⊤b. The matrix multiplying b is important it is called the pseudo-inverse of A,  enough to have a name and special notation: denoted by  A† := (A⊤A)−1A⊤. The solution via the normal equations using direct solvers amounts to the four steps laid out in the algorithm below. Algorithm: Least Squares via Normal Equations  (1) Form B = A⊤A and y = A⊤b. (2) Compute the Cholesky Factor, i.e., the lower triangular matrix G satisfying  B = GG⊤.  (3) Solve the lower triangular system Gz = y for z. (4) Solve the upper triangular system G⊤x = z for x.  The overall computational work for the first step of the algorithm is approximately mn2 floating point operations (flops); for step 2 it is n3/3 + O(n2) flops, while steps 3 and 4 cost O(n2) flops. This is another case where although operation counts are in general unsatisfactory for measuring true performance, they do deliver the essential result that the main cost here, especially when m ≫ n, is in forming the matrix B = A⊤A.  2. LEAST SQUARES AND THE NORMAL EQUATIONS  2.4. Data Fitting. Generally, data fitting problems arise as follows. We have observed data b and a model function that for any candidate model x provides predicted data. The task is to find x such that the predicted data match the observed data to the extent possible, by minimizing their difference in the least squares sense. In the linear case which we study here, the predicted data are given by Ax. In this context the assumption that A has full column rank does not impose a serious restriction: it just implies that there is no redundancy in the representation of the predicted data, so that for any vector ˆx ∈ Rn there is no other vector ˜x ∈ Rn such that A˜x = Aˆx.  Example 2.2. Consider fitting a given data set of m pairs (ti, bi) by a straight line.  Thus, we want to find the coefficients x1 and x2 of  v(t) = x1 + x2t,  such that v(ti) ≈ bi, i = 1, . . . , m. So n = 2 here, and  A =   t1 1 t2 1   ... ...   1 tm         .  The components of the normal equations are  B11 =  m (cid:88)  1 = m, B12 =  m (cid:88)  1 · ti =  m (cid:88)  ti,  i=1  i=1  i=1  B21 = B12, B22 =  m (cid:88)  ti · ti =  m (cid:88)  t2 i ,  y1 =  m (cid:88)  bi, y2 =  i=1 m (cid:88)  tibi.  i=1  i=1  i=1  This leads to a system of just two equations given by  mx1 +  (cid:32) m (cid:88)  ti  (cid:33)  x2 =  m (cid:88)  bi,  (cid:32) m (cid:88)  ti  (cid:33)  x1 +  i=1 (cid:32) m (cid:88)  t2 i  (cid:33)  x2 =  i=1 m (cid:88)  tibi.  i=1  i=1  i=1  The solution is written explicitly as the famous formula  i=1 bi − (cid:80)m i − ((cid:80)m i=1 t2 i=1 ti i − ((cid:80)m While in a typical course on statistics for the social sciences the regression formulas of Example 2.2 appear as if by magic, here they are a simple by-product of a general treatment – albeit, without the statistical significance.  (cid:80)m  (cid:80)m i=1 t2 i m (cid:80)m i=1 tibi − (cid:80)m m (cid:80)m i=1 t2  (cid:80)m  i=1 tibi i=1 ti)2 (cid:80)m i=1 bi i=1 ti)2  i=1 ti  x1 =  ,  m (cid:80)m  x2 =  ,  187  188  CHAPTER 6.  INTRODUCTION TO NUMERICAL OPTIMIZATION  2.5. Polynomial Data Fitting. Extending the linear regression formulas to a higher degree polynomial fit, v(t) ≡ pn−1(t) = x1 + x2t + . . . + xntn−1, is straightforward too. Writing for each data point v(ti) = x1 + x2ti + . . . + xntn−1 )x, the matrix A is the extension of the previously encountered Vandermonde matrix given by  i = (1, ti, . . . , tn−1  i   1 t1 1 t2   ... ...    1 tm−1  tm 1    t2 1 t2 2  tn−1 1 tn−1 2 ... tn−1 m−1 tn−1 m  · · · · · ... · · · · · ·         A =  t2 m−1 t2 m Note that the structure of the matrix A depends on our choice of the basis functions used to describe polynomials.  Here is a MATLAB function for best approximation by low order polynomials, using  the normal equations:  function coefs = lsfit (t , b , n ) % % function coefs = lsfit (t , b , n ) % % Construct coefficients of the polynomial of % degree at most n -1 that best fits data (t , b ) t = t (:) ; b = b (:) ; % make sure t and b are column vectors m = length ( t ) ; % long and skinny A A = ones (m , n ) ; for j =1: n -1  A (: , j +1) = A (: , j ) .* t ;  end % normal equations and solution B = A '* A ; y = A '* b ; coefs = B \ y ;  Example 2.3. Sample the function f (t) = cos(2πt) at 21 equidistant points on the interval [0, 1] and construct best fits by polynomials of degree at most n − 1 for each n = 1, 2, 3, 4, 5.  Following is an appropriate MATLAB script:  % data m = 21; tt = 0:1/( m -1) :1; bb = cos (2* pi * tt ) ; % find polynomial coefficients for n =1:5  coefs { n } = lsfit ( tt , bb , n ) ;  2. LEAST SQUARES AND THE NORMAL EQUATIONS  Figure 3. The first 5 best polynomial approximations to f (t) = cos(2πt) sampled at 0 : .05 : 1. The data values appear as red cir- cles. Clearly, p4 fits the data better than p2, which in turn is a better approximation than p0. Note p2j+1 = p2j.  end % Evaluate and plot t = 0:.01:1; z = ones (5 ,101) ; for n =1:5  z (n ,:) = z (n ,:) * coefs { n }( n ) ; for j =n -1: -1:1  z (n ,:) = z (n ,:) .* t + coefs { n }( j ) ;  end  end plot (t ,z , tt , bb , ' ro ') xlabel ( 't ') ylabel ( ' p_ {n -1} ')  The resulting approximants pn−1(t) are plotted in Figure 3. Note that here, due to symmetry, p1(t) = p0(t) and p3(t) = p2(t). So, the degree “at most n − 1” turns out to be equal to n − 2 rather than to n − 1 for odd values of n − 1.  189  190  CHAPTER 6.  INTRODUCTION TO NUMERICAL OPTIMIZATION  2.6. Data Fitting vs. Interpolation. Apparently, there is a seeming paradox hidden in our arguments, namely, that we attempt to minimize the residual for a fixed n, n < m, but refuse to simply increase n until n = m. Indeed, this would drive the residual to zero, the resulting scheme being a polynomial that interpolates the data. Why not just interpolate?!  The simple answer is that choosing n is part of our modeling efforts, and the ensuing least squares minimization problem is part of the solution process with n already fixed. But there are reasons for choosing n small in the first place. One is hinted at in Example 2.2, namely, that we are trying to find the trend in the data on a long time scale.  An additional reason for not interpolating the data values is that they may contain measurement errors. Moreover, we may want a model function that depends only on a few parameters xj for ease of manipulation, although we determine them based on all the given data. The solution of the least squares data fitting problem through solving the normal equations has the advantage of being straightforward and efficient. A linear least squares solver is implemented in the MATLAB backslash operator as well as in the MATLAB command polyfit. Replacing the last three lines in our function lsfit by the line  coefs = A \ b ;  would implement for the same purpose, albeit more enigmatically, an algorithm which in terms of roundoff error accumulation is at least as good. The routine polyfit is even easier to use, as it does not require forming the matrix and the right-hand-side vector; the input consists of the data points and the degree of the required polynomial.  2.7. Data Fitting in other norms. Before moving on let us also mention data  fitting in other norms.  (1) Using l1 we consider  min x  ∥b − Ax∥1.  Here we would be looking at finding coefficients x such that the sum of absolute values of the deviations is minimized. This norm is particularly useful if we need to automatically get rid of an undue influence of outliers in the data, which are data values that conspicuously deviate from the rest due to measurement error.  (2) Using l∞ we consider  min x  ∥b − Ax∥∞.  This is a min-max problem of finding the minimum over x1, x2, . . . , xn of the maximum data deviation. This norm is useful if the worst-case error in the approximation is important and must be kept in check.  Both l1 and l∞ best approximations lead to linear programming problems. They are significantly more complex than the problem that we are faced with using least squares yet are very important in many modern areas of application.  Note that if n = m, then we have one and the same solution regardless of the norm used, because with all three norms the minimum is obtained with x = A−1b, which yields a zero residual.  2. PROBLEMS  But when n < m these different norms usually yield significantly different best  approximations.  The least squares approximation (which is the approximation associated with the l2- norm) is not only the simplest to calculate, it also has a variety of beautiful mathematical properties.  Problems  (1) The following data were obtained in connection with the values of a certain  fictitious material property:  t 0.0 1.0 0.5 b 0.9 1.01 1.05 0.97 0.98 0.95 0.01 −0.1 0.02 −0.1 0.0  0.1  0.2  0.3  0.4  0.6  0.7  0.8  0.9  It was then hypothesized that the underlying material property is a piecewise constant function with one break point (i.e., two constant pieces). (a) Plot the data points and decide where (approximately) the break point  should be.  (b) Find the piecewise constant function which best fits the data by least  squares, and plot it, too.  (2)  (a) Evaluate the function  f (t) = .05 sin(1000t) + .5 cos(πt) − .4 sin(10t)  at the 101 points given by 0 : .01 : 1. Plot the resulting broken line interpolant.  (b) In order to study the slow scale trend of this function, we wish to find a low degree polynomial (degree at most 6) that best approximates f in the least squares norm at the above 101 data points. By studying the figure from part (a) find out the smallest n that would offer a good fit in this sense. (Try to do this without further computing.)  (c) Find the best approximating polynomial v of degree n and plot it together  with f . What are your observations?  (3) Let us synthesize data in the following way. We start with the cubic polynomial  q(t) = −11 +  55 3  t −  17 2  t2 +  7 6  t3.  Thus, n = 4. This is sampled at 33 equidistant points between 0.9 and 4.1. Then we add to these values 30% noise using the random number generator randn in MATLAB, to obtain “the data” which the approximations that you will construct “see”. From here on, no knowledge of q(t), its properties, or its values anywhere is assumed: we pretend we don’t know where the data came from!  Your programming task is to pass through these data three approximations: (a) An interpolating polynomial of degree 32. This can be done using the MATLAB function polyfit. You don’t need to know how such interpo- lation works for this exercise.  (b) An interpolating cubic spline using the MATLAB function spline. (c) A cubic polynomial which best fits the data in the l2 sense, obtained by  our function lsfit.  191  192  CHAPTER 6.  INTRODUCTION TO NUMERICAL OPTIMIZATION  Plot the data and the obtained approximations. Which of these approximations make sense? Discuss.  (4) Often in practice, an approximation of the form  u(t) = γ1eγ2t  is sought for a data fitting problem, where γ1 and γ2 are constants. Assume given data (t1, z1), (t2, z2), . . ., (tm, zm), where zi > 0, i = 1, 2, . . . , m, and m > 0. (a) Explain in one brief sentence why the techniques introduced in the present  chapter cannot be directly applied to find this u(t).  (b) Considering instead  v(t) = ln u(t) = (ln γ1) + γ2t,  it makes sense to define bi = ln zi, i = 1, 2, . . . , m, and then find coefficients x1 and x2 such that v(t) = x1 + x2t is the best least squares fit for the data (t1, b1), (t2, b2), . . ., (tm, bm). Using this method, find u(t) for the data  i ti zi  1 0.0 e0.1  2 1.0 e0.9  3 2.0 e2  Bibliography  [1] E.W. Cheney, D. Kincaid, Numerical Mathematics and Computing. Sixth Edition. Brooks/Cole,  2008.  [2] U. M. Ascher, C. Greif, A First Course in Numerical Methods. SIAM Society for Industrial and  Applied Mathematics, Philadelphia.  [3] J. H. Mathews, K. D. Fink, Numerical Methods Using MATLAB. Pearson Prentice Hall.  193

JSON: