Precisión finita, redondeo y condicionamiento

Cálculo con números de precisión finita

En los capítulos anteriores tratamos los algoritmos numéricos como transformaciones definidas sobre estados, parámetros y pasos temporales. Esa descripción bastó para estudiar la discretización, la estabilidad, las tolerancias y la rigidez. Falta incorporar una limitación más elemental. Un computador ordinario no representa todos los números reales y, por tanto, tampoco ejecuta exactamente todas las operaciones de la aritmética real.

Eso no significa que los computadores sean imprecisos en un sentido vago. La aritmética de punto flotante está definida mediante formatos y reglas precisas. El estándar IEEE 754 especifica las representaciones y los modos de redondeo, además del comportamiento ante situaciones excepcionales (IEEE Computer Society, 2019). En un entorno típico de cálculo científico, el tipo float usa el formato binario de 64 bits, llamado binary64.

Un resultado numérico puede diferir del resultado matemático por varias razones. Los datos pueden ser aproximados y el modelo puede ser una idealización. El método puede discretizar una ecuación continua y la propia aritmética puede redondear las operaciones intermedias. Estas fuentes de discrepancia no deben agruparse bajo la palabra «error» como si fueran equivalentes.

El capítulo parte del error que introduce representar números finitos y operar con ellos. Sobre esa base distingue la sensibilidad de la respuesta correcta frente a pequeñas perturbaciones de los datos del error adicional que introduce el algoritmo concreto con que se resuelve el problema.

La distinción clásica entre condicionamiento del problema y estabilidad del algoritmo responde a esas dos cuestiones. Un problema puede ser intrínsecamente sensible aunque el algoritmo sea excelente. También puede ocurrir lo contrario, cuando el problema resiste bien las perturbaciones pequeñas y una formulación numérica concreta destruye información útil. Higham (2002) organiza buena parte del análisis numérico moderno en torno a la separación entre la sensibilidad del problema y los errores que introduce el cálculo.

Qué significa representar un número en punto flotante

Un número binario normalizado de punto flotante puede escribirse conceptualmente como $\pm m\,2^e$, donde $m$ contiene un número finito de bits significativos y $e$ determina la escala. En binary64 se almacenan un bit de signo, once bits de exponente y cincuenta y dos bits explícitos de fracción. Los números normales tienen además un bit significativo implícito, de modo que la precisión efectiva es de 53 bits (Goldberg, 1991; IEEE Computer Society, 2019).

Eso no equivale a «quince decimales exactos» en todas las circunstancias. La precisión relativa se mide en la base binaria y la separación absoluta entre números representables crece con su magnitud. Cerca de 1 esa separación es de $2{,}2\times10^{-16}$, cerca de 1000 es de $1{,}1\times10^{-13}$ y cerca de $10^{16}$ vale 2.

Dos constantes deben distinguirse, porque con frecuencia reciben el mismo nombre en distintos textos. En binary64, la distancia entre 1 y el siguiente número representable mayor que 1 es

Con redondeo al número representable más próximo, el redondeo unitario suele definirse como

En condiciones ordinarias, sin desbordamiento, subdesbordamiento problemático ni otros casos especiales, una operación elemental puede modelarse como

donde $\circ$ designa una operación básica y $\operatorname{fl}$ el resultado redondeado. La expresión da una cota local útil para el análisis, aunque no afirma que todos los errores tengan exactamente la magnitud $u$. Goldberg (1991) explica en detalle por qué estos pequeños errores pueden combinarse de maneras muy distintas según cómo se formule el cálculo. El cuadro 27.1 resume las propiedades de binary64 que usaremos.

Cuadro 27.1. Propiedades del formato binary64.

Los números subnormales extienden la representación por debajo del menor número normal, aunque ya no conservan la misma estructura relativa de espaciado. Por eso el modelo multiplicativo anterior no debe aplicarse mecánicamente en los extremos del dominio.

El formato finito introduce además límites de rango. Si una magnitud excede el máximo representable, el resultado puede convertirse en infinito según las reglas del estándar. En el extremo opuesto, los subnormales permiten acercarse gradualmente a cero, con menor precisión relativa. Por eso hay que separar el desbordamiento y el subdesbordamiento del redondeo ordinario. Todos son consecuencias de la representación finita, pero generan patrones de error distintos y exigen diagnósticos diferentes. Un algoritmo puede ser estable respecto del redondeo ordinario y fallar, sin embargo, porque una variable intermedia sale del rango representable.

El redondeo altera algunas propiedades algebraicas

En la aritmética real, la suma es asociativa, es decir, $(a+b)+c=a+(b+c)$. En punto flotante los resultados intermedios se redondean. Consideremos $a=10^{16}$, $b=-10^{16}$ y $c=1$. En binary64, $(a+b)+c=1$, porque $a+b$ se representa exactamente como cero y luego se suma uno. En cambio, $a+(b+c)=0$. La suma $b+c$ se redondea a $-10^{16}$, porque a esa escala la separación entre números representables es 2. La información que aportaba $c$ desaparece antes de la suma final.

El ejemplo muestra que las operaciones se realizan en un conjunto finito y que cada resultado se proyecta de nuevo sobre ese conjunto. Las reglas de IEEE 754 están bien definidas, aunque la estructura resultante no conserva todas las identidades de los números reales.

Por eso, reordenar una suma larga o paralelizarla puede modificar los últimos dígitos del resultado. Si una conclusión científica depende de que esas diferencias no existan, afirmar que el código es determinista no resuelve el problema. Hay que evaluar la sensibilidad de la cantidad calculada y del algoritmo usado.

La sumación compensada de Kahan (1965) es un ejemplo clásico de algoritmo diseñado para conservar parte de la información que una suma secuencial ingenua pierde. Si a 1 se le suman un millón de términos iguales a $10^{-16}$ en ese orden, la suma ingenua devuelve exactamente 1, porque cada término queda por debajo de la mitad de la separación entre números cercanos a 1. La suma compensada devuelve $1{,}0000000001$, el valor correcto. El método no hace exacta cualquier suma, pero muestra un principio que reaparecerá a lo largo del libro. Dos algoritmos algebraicamente equivalentes pueden comportarse de manera muy distinta en precisión finita.

Cancelación al restar números cercanos

Consideremos la función $g(x)=\sqrt{x+1}-\sqrt{x}$, que puede racionalizarse como

Las dos expresiones son exactamente equivalentes sobre los reales positivos, aunque numéricamente se comportan de modo distinto. Para $x=10^{16}$, la forma directa evaluada en binary64 da 0. El número $10^{16}+1$ no es representable y se redondea a $10^{16}$, de modo que las dos raíces coinciden. La forma racionalizada da $g(10^{16})\approx5\times10^{-9}$, que coincide con una evaluación de alta precisión en todos los dígitos mostrados.

Cuadro 27.2. Evaluación de $g(x)=\sqrt{x+1}-\sqrt{x}$ en binary64 con la forma directa y con la forma racionalizada. El error relativo se calculó frente a una evaluación con 50 dígitos.

La figura 27.1 muestra el error relativo de ambas formas en todo el intervalo. El de la forma directa crece aproximadamente en proporción a $x$ hasta llegar a 1, mientras que el de la forma racionalizada permanece en el nivel del redondeo unitario.

El diagnóstico debe formularse con cuidado. El problema matemático no se vuelve mal condicionado por usar una mala expresión para evaluarlo. Para valores grandes de $x$, su número de condición relativo respecto de $x$ se acerca a $1/2$, de modo que la función responde con moderación a perturbaciones relativas de su argumento. La pérdida de precisión procede de la formulación algorítmica. El problema es poco sensible. Este es el primer ejemplo nítido de la distinción que buscamos, porque un problema razonablemente condicionado se resuelve aquí con un procedimiento numéricamente inadecuado.

Condicionamiento y sensibilidad de la respuesta correcta

Sea una función escalar $y=f(x)$. Para una perturbación pequeña $\Delta x$, la aproximación lineal da $\Delta y\approx f'(x)\Delta x$. Una medida local de sensibilidad absoluta es $|f'(x)|$. Cuando tiene sentido comparar cambios relativos, se define

Este número indica, aproximadamente, cuántas unidades de error relativo en la salida pueden producirse por cada unidad de perturbación relativa en la entrada, para perturbaciones suficientemente pequeñas.

La definición debe interpretarse con cuidado. Depende de la variable elegida y de la noción de perturbación, además de la norma en los problemas vectoriales. El número de condición describe la pregunta tal como fue formulada y puede cambiar si la pregunta se formula de otro modo.

Para un sistema lineal $Ax=b$ con $A$ invertible, el número de condición asociado a una norma subordinada es

Si $\kappa(A)$ es cercano a uno, pequeñas perturbaciones relativas de los datos producen cambios relativos comparables en la solución. Si es muy grande, existe alguna dirección de perturbación que puede amplificarse notablemente. El número de condición no predice que toda perturbación se amplifique exactamente por $\kappa(A)$. Mide una sensibilidad potencial y entra en cotas que requieren hipótesis adicionales. Higham (2002) desarrolla versiones por norma y por componentes, porque una única cifra puede ocultar estructuras importantes de escala y de signo.

Para ver de dónde procede el condicionamiento, supongamos que $A$ se mantiene fija y que solo se perturba $b$. Como $A(x+\Delta x)=b+\Delta b$, al restar $Ax=b$ se obtiene $A\Delta x=\Delta b$, es decir, $\Delta x=A^{-1}\Delta b$. Por la submultiplicatividad de la norma, $\|\Delta x\|\le \|A^{-1}\|\,\|\Delta b\|$. De $b=Ax$ se obtiene además $\|b\|\le\|A\|\,\|x\|$. Al combinar ambas relaciones resulta

La desigualdad es una cota y no predice el error exacto para cada perturbación. Su utilidad es mostrar por qué un $\kappa(A)$ grande permite amplificar variaciones pequeñas de los datos. Cuando también se perturba $A$, la expresión requiere términos adicionales y condiciones que eviten acercarse demasiado a la singularidad (Higham, 2002).

Un sistema casi dependiente

Consideremos

Para cualquier $\delta\ne0$ la solución exacta es $x=(1,1)^{\mathsf T}$. Cuando $\delta$ es pequeño, las dos ecuaciones corresponden a rectas casi paralelas. Si solo se modifica el segundo dato de $b$ en una cantidad $\delta$, con $\widetilde b=(2,\ 2+2\delta)^{\mathsf T}$, la nueva solución es $\widetilde x=(0,2)^{\mathsf T}$. La figura 27.2 usa $\delta=0{,}02$.

Tomemos ahora $\delta=10^{-8}$. El número de condición en norma 2 es de unos $4{,}0\times10^8$ y la perturbación relativa aplicada a $b$ es apenas

mientras que el cambio relativo de la solución es

Una modificación de unas pocas partes por mil millones en los datos duplica una componente de la solución y elimina la otra. El efecto ocurre en la aritmética real exacta, sin intervención del redondeo, porque el problema está mal condicionado. Para esta familia, $\kappa_2(A_\delta)$ crece aproximadamente como $4/\delta$ cuando $\delta$ es pequeño (figura 27.3).

Error hacia adelante, error hacia atrás y residuo

Supongamos que el problema verdadero es $Ax=b$ y que un algoritmo devuelve una aproximación $\widehat x$. El error hacia adelante compara directamente la respuesta calculada con la verdadera, mediante $\|\widehat x-x\|$. El residuo, $r=b-A\widehat x$, mide cuánto incumple $\widehat x$ la ecuación original. El error hacia atrás plantea otra pregunta, la de cuánto habría que perturbar los datos para que $\widehat x$ fuera la solución exacta de un problema cercano. Wilkinson (1963) convirtió esta inversión de perspectiva en una herramienta central del análisis numérico.

Volvamos al ejemplo con $\delta=10^{-8}$ y tomemos $\widehat x=(0,2)^{\mathsf T}$. Respecto del sistema original, el cambio relativo de la solución es de orden uno. El residuo relativo, en cambio, es de unos $3{,}54\times10^{-9}$. Es decir, $\widehat x$ resuelve exactamente un sistema cuyo vector $b$ difiere del original apenas en esa escala relativa. Un algoritmo perfectamente estable hacia atrás podría devolver una solución con gran error hacia adelante cuando el problema está mal condicionado.

El resultado corrige una interpretación frecuente, según la cual un residuo pequeño garantiza una solución precisa. Un residuo pequeño indica que la solución satisface muy bien una ecuación cercana. Para convertir esa afirmación en una cota del error hacia adelante hay que incorporar el condicionamiento (Higham, 2002; Turing, 1948). En términos esquemáticos, cuando se cumplen las condiciones de la teoría perturbacional correspondiente,

Problema sensible y algoritmo inestable

Los ejemplos pueden ordenarse ahora. En la diferencia $\sqrt{x+1}-\sqrt{x}$, el problema matemático está razonablemente condicionado para $x$ grande, aunque la evaluación directa sufre cancelación. La forma racionalizada obtiene mucha mayor precisión, así que el diagnóstico principal es algorítmico. En el sistema casi dependiente, en cambio, una perturbación diminuta de los datos puede modificar sustancialmente la solución incluso con aritmética exacta. Ningún algoritmo puede recuperar en la solución una precisión que los datos, tal como están planteados, no contienen de manera estable. El diagnóstico principal es el condicionamiento.

Ambos problemas pueden coexistir. Un sistema mal condicionado puede resolverse además con un algoritmo inestable, que amplifica todavía más los errores. También puede resolverse con un método estable y conservar una solución sensible. Por eso usar más precisión no es una respuesta universal. Más bits reducen el error de redondeo, pero no modifican el condicionamiento matemático del problema original.

La elección de la formulación también importa. En los mínimos cuadrados, resolver las ecuaciones normales eleva aproximadamente al cuadrado el número de condición de la matriz de diseño. Las factorizaciones QR o SVD no cambian los datos y evitan esa fuente adicional de sensibilidad numérica. La estructura algebraica del algoritmo puede, así, conservar o degradar la información disponible.

Formulación, escala y número de condición

El número de condición depende de la formulación. En un sistema lineal, $\kappa(A)$ depende de la norma usada y puede cambiar si se reescalan las variables. Supongamos que una componente del estado se expresa en segundos y otra en milisegundos. Multiplicar una coordenada por mil modifica los coeficientes de la matriz y puede cambiar mucho su número de condición, aunque el fenómeno sea el mismo. La matriz identidad de orden dos, cuyo número de condición es 1, pasa a tener número de condición 1000 si una de sus variables se mide en una unidad mil veces menor.

La observación tiene dos consecuencias. Antes de declarar que un modelo está mal condicionado hay que examinar si la dificultad procede de unidades o escalas muy desiguales que pueden normalizarse de manera científicamente legítima. Reducir el número de condición mediante un cambio de variables, por otra parte, no garantiza que se haya eliminado una sensibilidad real, porque puede tratarse solo de una métrica más adecuada para comparar perturbaciones.

En álgebra lineal numérica, el precondicionamiento transforma el problema en otro equivalente que resulte más favorable para un algoritmo iterativo. Para resolver $Ax=b$ puede introducirse, por ejemplo, una matriz invertible $M$ y resolverse $M^{-1}Ax=M^{-1}b$. La solución $x$ no cambia, aunque la geometría numérica del problema sí puede cambiar. Un buen precondicionador aproxima la estructura difícil de $A$ con un costo muy inferior al de resolver exactamente el problema original. Más adelante aparecerá una idea semejante en los sistemas dinámicos, donde transformar una representación puede facilitar el cálculo sin alterar el objeto matemático estudiado.

También hay que distinguir el análisis por norma del análisis por componentes. Una norma resume las perturbaciones en una sola magnitud, lo que es útil, pero puede ocultar que algunas componentes se conocen con precisiones muy distintas o que ciertos ceros tienen significado estructural. Por eso dos problemas con el mismo número de condición por norma pueden responder de forma distinta a las perturbaciones físicamente plausibles. El condicionamiento debe leerse siempre junto con la clase de perturbaciones que se consideran admisibles.

Precisión finita en sistemas dinámicos

Los sistemas dinámicos agregan una dificultad, porque el resultado de una operación se convierte en entrada del paso siguiente. Un error pequeño puede amortiguarse o amplificarse según la dinámica y el algoritmo. En un sistema contractivo, las trayectorias cercanas tienden a aproximarse y el efecto de las perturbaciones pequeñas puede disminuir. En un sistema con direcciones expansivas, una perturbación diminuta puede crecer. Más adelante formalizaremos esa expansión mediante los exponentes de Lyapunov.

Hay que evitar una inferencia prematura. Que dos simulaciones diverjan al cambiar los últimos bits de la condición inicial no demuestra caos. La divergencia puede deberse a una inestabilidad numérica o a una propiedad real de la dinámica, entre otras causas. Cada posibilidad requiere controles distintos.

La precisión finita tiene además una consecuencia lógica. Si una iteración determinista opera sobre un conjunto finito de estados representables y permanece dentro de él, tarde o temprano visita un estado ya recorrido y desde ese momento la secuencia es periódica. Una trayectoria digital infinita no puede, por tanto, reproducir literalmente para siempre una órbita matemática aperiódica. La iteración logística $x_{n+1}=4x_n(1-x_n)$ lo muestra con claridad. Desde $x_0=0{,}3$, con el tipo float16 de NumPy la secuencia entra en un ciclo de período 40 después de 35 pasos. Con float32, el ciclo tiene período 4344 y aparece después de 224 pasos.

Esto no vuelve inútil la simulación del caos. Obliga a definir qué propiedades se espera reproducir, como distribuciones, exponentes, tiempos de predicción o estadísticas invariantes. La coincidencia puntual indefinida con una trayectoria ideal queda fuera de alcance. También obliga a repetir los cálculos con distintas precisiones, algoritmos y pasos cuando el resultado científico podría depender de esos detalles.

Controles reproducibles para este capítulo

Los ejemplos anteriores pueden reproducirse con pocas líneas de código. El propósito es convertir las afirmaciones numéricas del capítulo en objetos verificables.

import math
import numpy as np

f = np.finfo(float)
assert f.eps == 2.0**-52          # distancia entre 1 y el siguiente float
u = f.eps / 2                     # redondeo unitario
assert (1e16 - 1e16) + 1.0 == 1.0
assert 1e16 + (-1e16 + 1.0) == 0.0

x = 1e16
directa = math.sqrt(x + 1) - math.sqrt(x)
racionalizada = 1 / (math.sqrt(x + 1) + math.sqrt(x))
print(directa, racionalizada)

Otro bloque cuantifica el condicionamiento.

import numpy as np

delta = 1e-8
A = np.array([[1.0, 1.0],
              [1.0, 1.0 + delta]])
b = np.array([2.0, 2.0 + delta])
b2 = b + np.array([0.0, delta])

x1 = np.linalg.solve(A, b)
x2 = np.linalg.solve(A, b2)

print(f"{np.linalg.cond(A, 2):.3e}")
print(x1, x2)
print(f"{np.linalg.norm(b2 - b) / np.linalg.norm(b):.3e}")
print(f"{np.linalg.norm(x2 - x1) / np.linalg.norm(x1):.3f}")

El bloque inicial comprueba con aserciones el valor de $\varepsilon$ y la pérdida de asociatividad. Luego imprime 0.0 5e-09, la cancelación completa de la forma directa junto al valor correcto de la racionalizada. El bloque del sistema casi dependiente imprime el número de condición $4{,}000\times10^8$ y las dos soluciones $(1,1)$ y $(0,2)$. Imprime también la perturbación relativa de los datos, que vale $3{,}536\times10^{-9}$, junto al cambio relativo de la solución, igual a 1.

El interés de estas pruebas no está en que Python «confirme» la teoría. La teoría permite anticipar qué debe ocurrir y la ejecución comprueba que el código y las cifras del texto son consistentes con esa derivación.

Consolidación

Este capítulo añade una capa al esquema de error construido desde el capítulo 6. Un resultado numérico depende del modelo elegido y del método que lo aproxima, así como de los datos y de la aritmética concreta en que se ejecuta el cálculo. Estas capas pueden interactuar, aunque no deben confundirse.

La aritmética de punto flotante trabaja con un subconjunto finito de números y redondea las operaciones de acuerdo con reglas definidas. Cada redondeo elemental puede ser pequeño y, aun así, una secuencia de operaciones mal formulada puede destruir información mediante cancelación o acumulación desfavorable. Por eso operaciones algebraicamente equivalentes pueden tener comportamientos numéricos muy distintos.

El condicionamiento pertenece al problema formulado, porque mide cuánto puede cambiar la respuesta exacta cuando cambian ligeramente sus datos. La estabilidad pertenece al procedimiento, porque evalúa si el algoritmo introduce errores mucho mayores que los necesarios para explicar su resultado como solución de un problema cercano. La cancelación puede destruir cifras en un problema bien condicionado, mientras que un problema mal condicionado puede producir un gran error hacia adelante con un residuo y un error hacia atrás pequeños. Aumentar la precisión reduce el redondeo, pero no corrige un problema intrínsecamente mal condicionado. El cuadro 27.3 resume la pregunta que responde cada concepto.

Cuadro 27.3. Pregunta principal que responde cada concepto del capítulo.

En los sistemas dinámicos, los errores de representación se propagan según la geometría de la evolución. Su crecimiento, considerado aisladamente, no demuestra caos. Estas distinciones serán indispensables en los barridos paramétricos. Una transición aparente puede ser una propiedad del sistema, un efecto de la resolución de la cuadrícula, una inestabilidad del integrador o la amplificación de pequeñas perturbaciones en una región mal condicionada. El capítulo 28 desarrollará procedimientos para explorar parámetros sin convertir artefactos numéricos en cambios de régimen del modelo.

Ejercicios

Ejercicio 27.1. (a) Calcule la separación entre números binary64 consecutivos cerca de 1, de 1000 y de $10^{16}$. (b) Explique por qué la comparación 1e16 + 1 == 1e16 es verdadera en Python. (c) ¿Cuál es el mayor entero $n$ tal que todos los enteros entre 0 y $n$ son exactamente representables en binary64?

Ejercicio 27.2. A 1 se le suman, en ese orden, un millón de términos iguales a $10^{-16}$. (a) ¿Qué devuelve la suma secuencial ingenua y por qué? (b) ¿Qué ocurre si los términos pequeños se suman primero? (c) Implemente la suma compensada de Kahan y compare su resultado con el de math.fsum.

Ejercicio 27.3. La ecuación $x^2-10^8x+1=0$ tiene raíces cercanas a $10^8$ y a $10^{-8}$. (a) Calcule la raíz pequeña en binary64 con la fórmula habitual $x=\left(-b-\sqrt{b^2-4ac}\right)/(2a)$. (b) Calcúlela como $c/(a\,x_{\mathrm{grande}})$. (c) ¿Qué operación produce el error y qué nombre recibe el fenómeno?

Ejercicio 27.4. Calcule el número de condición relativo $\kappa_f(x)$ de (a) $f(x)=\sqrt{x+1}-\sqrt{x}$ para $x$ grande, (b) $f(x)=x-1$ cerca de $x=1$ y (c) $f(x)=e^x$. ¿En cuál de los tres casos la dificultad pertenece al problema mismo?

Ejercicio 27.5. Para la matriz $A_\delta$ del capítulo, (a) calcule $\kappa_2(A_\delta)$ con $\delta$ igual a $10^{-2}$, $10^{-4}$ y $10^{-8}$ y compárelo con $4/\delta$. (b) Si el vector $b$ se conoce con un error relativo de $10^{-9}$ y $\delta=10^{-8}$, ¿qué error relativo en la solución admite la cota del condicionamiento? (c) Explique por qué resolver el sistema en cuádruple precisión no mejora la respuesta de (b).

Ejercicio 27.6. Itere la función logística $x_{n+1}=4x_n(1-x_n)$ desde $x_0=0{,}3$ con los tipos float16 y float32 de NumPy. (a) Detecte el primer estado repetido y calcule la longitud del transitorio y el período del ciclo. (b) ¿Por qué la repetición es inevitable? (c) ¿Qué propiedades de una dinámica caótica puede aspirar a reproducir una simulación en precisión finita?

def suma_kahan(valores):
    s = 0.0
    c = 0.0
    for a in valores:
        y = a - c
        t = s + y
        c = (t - s) - y
        s = t
    return s

Referencias

Goldberg, D. (1991). What every computer scientist should know about floating-point arithmetic. ACM Computing Surveys, 23(1), 5–48. https://doi.org/10.1145/103162.103163

Higham, N. J. (2002). Accuracy and stability of numerical algorithms (2.ª ed.). Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898718027

IEEE Computer Society. (2019). IEEE standard for floating-point arithmetic (IEEE Std 754-2019). https://doi.org/10.1109/IEEESTD.2019.8766229

Kahan, W. (1965). Further remarks on reducing truncation errors. Communications of the ACM, 8(1), 40. https://doi.org/10.1145/363707.363723

Turing, A. M. (1948). Rounding-off errors in matrix processes. The Quarterly Journal of Mechanics and Applied Mathematics, 1(1), 287–308. https://doi.org/10.1093/qjmam/1.1.287

Wilkinson, J. H. (1963). Rounding errors in algebraic processes. Her Majesty's Stationery Office.