In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
В этой лабораторной мы попробуем оценить точность результата вычисления по внешне невинной формуле в арифметике с плавающей запятой. Обычно делая арифметические операции мы не задумываемся над точностью результата, однако даже одна операция может значительно ухудшить точность. Рассмотрим несколько примеров.
Чтобы оценить точность вычислений, нам необходимо знать точный ответ. Так как мы только хотим проиллюстрировать одну из проблем с потерей точности, то мы можем сами сконструировать формулу так, чтобы точных ответ был заранее известен. Мы будем вычислять значение функции, которая эквивалента в точной арифметике тождественному преобразованию, т.е. точный результат вычисления фукнции совпадает с аргументом функции: $$f(x)=x.$$ Однако в арифметике с плавающей запятой не все свойства арифметических операций верны, поэтому результат будет получаться с ошибкой. Мы будем вычислять относительную ошибку по обыкновенной формуле: $$\sigma(x_0,x)=\frac{|x_0-x|}{|x_0|},$$ где $x_0$ - точное значение, а $x$ - приближенное.
In [2]:
def relative_error(x0, x): return np.abs(x0-x)/np.abs(x0)
Исходные данные нам даны с точностью представления вещественных чисел числами с плавающей запятой, т.е. с половиной машинной точности. Мы будем проводить вычисления в арифметике двойной точности.
In [3]:
eps=np.finfo(np.double).eps
print("Машинная точность:",eps)
Рассмотрим функцию, выполняющую $n$ делений на число $d$, а затем $n$ умножений на $d$. Ясно, что в точной арифметике эта функция не должна изменять аргумент.
In [4]:
def f_div_mult(x, d=np.pi, n=52):
for k in range(n): x=x/d
for k in range(n): x=x*d
return x
Создадим набор тестовых значений и проверим точность вычислений на этом наборе.
In [5]:
x0=np.logspace(-4,4,100,dtype=np.double)
x=f_div_mult(x0)
err=relative_error(x0, x)
print("Ошибки",err[:4],"...")
Гораздо лучше результаты воспринимаются на картинке. При оценки точности мы всегда выводим относительную ошибку в логарифмической шкале.
In [6]:
def plot_error(x0,err):
mask=np.logical_and(err>0,err<np.inf)
plt.loglog(x0[mask],err[mask],".k")
plt.loglog(x0,[eps]*len(err),"--r") # машинная точность для сравнения
plt.xlabel("$Значение\;аргумента$")
plt.ylabel("$Относительная\;погрешность$")
plt.show()
plot_error(x0,err)
Мы видим, что погрешность в несколько раз больше машинной точности, обозначенной красной пунктирной линией, однако погрешность остается пренебрежимо малой. Рассмотрим другую функцию, очень похожу на предыдущую, однако теперь мы будем сначала извлекать корни, а затем возводить результат в квадрат.
In [7]:
def f_sqrt_sqr(x, n=52):
for k in range(n): x=np.sqrt(x)
for k in range(n): x=x*x
return x
Проведем вычисления и построим график ошибки на том же наборе значений.
In [8]:
x=f_sqrt_sqr(x0)
err=relative_error(x0, x)
plot_error(x0, err)
Ошибка почти 100%! Как такое могло произойти? Не возникло исключений, нет сообщений об ошибках, но результат не содержит ни одной верной цифры.
Большинство людей предполагают, что причина столь большой ошибки заключается в том, что квадратный корень нельзя вычислить точно. Чтобы проверить это утверждение мы снова вычислим 52 квадратных корня и 52 квадрата, но изменим порядок вычислений на следующий: корень, квадрат, корень, квадрат,..
In [9]:
def f_sqrt_sqr_interleave(x, n=52):
for k in range(n):
x=np.sqrt(x)
x=x*x
return x
x=f_sqrt_sqr_interleave(x0)
err=relative_error(x0, x)
plot_error(x0, err)
Теперь мы получаем результат примерно с машинной точностью. Как можно обьяснить тот факт, что выполняя одни и те же операции в разном порядке, можно получить на много порядков отличающиеся ошибки?
Рассмотрим приближенное значение $x$ некоторой величины $x_0$ с относительной погрешностью $\delta$, тогда $x$ выражается через точное значение и погрешность следующим образом: $$x=x_0(1+\delta).$$ Рассмотрим теперь некоторую степень $\alpha$ значения $x$: $$x^\alpha=x_0^\alpha(1+\delta)^\alpha\approx x_0^\alpha(1+\alpha\delta),$$ где мы воспользовались эквивалентностью $(1+\delta)^\alpha\sim 1+\alpha\delta$ для бесконечно малых $\delta$. Второе слагаемое в скобках равно погрешностьи вычисления степени, т.е. при вычислении степени $\alpha$ погрешность аргумента увеличивается в $\alpha$ раз. Следовательно при вычислении квадратов относительная погрешность быстро увеличивается, а при вычислении квадратных корней - уменьшается.
Если при вычислении корней получен ответ с машинной точностью $\epsilon$, то в результате вычисления 52 квадратов погрешность достигнет значения $2^{52}\epsilon$, что неприемлемо. Следовательно нужно вычислить квадртаные корни с точностью превосходящей $\epsilon$. Однако сохранить число с плавающей запятой точнее чем с половиной машинной точности нельзя, следовательно необходимо изменить представления чисел. Создание новой реализации чисел с плавающей запятой весьма трудоемко, поэтому создадим свою реализацию на основе имеющихся чисел с плавающей запятой, но будем хранить не само число, а некоторую функцию числа. Например, заметим, что при вычислении квадратного корня результат быстро стремится к 1, поэтому мы можем представлять числа в виде $x=1+\zeta$, где $\zeta$ мало. Так как числа с плавающей запятой могут выражать очень малые величины, много меньше машинной точности, то такое представление позволит повысить точность вычислений. При сложении $\zeta$ с 1 результат будет округлен до 1, поэтому мы нигода не будем явно выполнять это сложение, а будем хранить только $\zeta$. Однако все арифметические операции должны выполняться так, как если бы они выполнялись над $x$, а не над $\zeta$. Например, так как результат перемножения $x=1+\zeta$ и $y=1+\eta$ равен $x\cdot y=1+\zeta+\eta+\zeta\cdot\eta$, то умножение на нашем представлении чисел определено не как обычное умножение, а как $$(\zeta,\eta)\mapsto \zeta+\eta+\zeta\cdot\eta.$$
In [10]:
# Пример собственной реализации чисел на основе чисел с плавающей запятой
class MyNumber(object):
def __init__(self, zeta):
"""Конструктор принимает zeta, но обьект соответствует числу x=1+zeta."""
self.zeta=zeta
def __str__(self):
"""На экран выводится значение x, которое может быть менее точно,
чем храниемое значение."""
return "{}".format(self.to_float())
def from_float(x):
"""Создает число со значением, равным x."""
return MyNumber(x-1)
def to_float(self):
"""Преобразует число в формат с плавающей запятой"""
return self.zeta+1
def __mul__(self, other):
"""Перезагрузка операции умножения."""
return MyNumber(self.zeta+other.zeta+self.zeta*other.zeta)
#def __mul__(self, other):
# """Так делать нельзя. Как сделать правильно?"""
# return MyNumber.from_float(self.to_float()*other.to_float())
In [11]:
print("Число с плавающей запятой:", np.pi)
pi=MyNumber.from_float(np.pi)
print("Наше представление числа: ", pi)
print("Квадрат в арифметике с плавающей запятой:", np.pi*np.pi)
print("Квадрат в нашем представлении: ", pi*pi)
Подумайте, как еще можно было представить числа? Какое преобразование переводи возведение в степень в умножение? Вместо хранения самого числа $x$ можно хранить $y$, такое что $x=\ln y$, тогда мы фактически получаем новое представление для чисел $x$, однако нам не нужно заново переписывать библиотеки для работы с числами, так как $y$ храниться в виде обычного числа с плавающей запятой. Убедитесь, что в таком представлении при вычислении степеней предельная погрешность почти не изменяется. Если это представление решает проблемы с погрешностью, почему такое представление не используется массово? В исходной формулировке задачи на вход поступает число $x$ в обычном виде числа с плавающей запятой, в таком же виде мы хотим иметь число на выходе вычислений. Какова погрешность преобразования чисел между форматами? Получаем ли мы выигрыш в точности с учетом преобразования форматов?
В машинном обучение часто используется следующий вариант LogSumExp гладкого максимума:
$$LSE(x) = \log(1+e^x).$$В данном случае $LSE(x)\approx \max(0,x)$, однако функция дифференциируема, что важно для использования этой функции в искуственных нейронных сетях. График функции действительно напоминает максимум, но сглажен.
In [12]:
def LSE(x):
return np.log(1+np.exp(x))
def max0(x):
return np.maximum(x, 0)
xs = np.linspace(-10,10,1000)
plt.plot(xs, max0(xs))
plt.plot(xs, LSE(xs))
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.legend(["$\max(0,x)$", "$LSE(x)$"])
plt.show()
Диапазон изменения аргумента однако может быть значительно шире представленного на графике. Даже для умеренных значениях аргумента наш способ вычисления функции $LSE$ в числах с плавающей запятой не приводит к успеху.
In [13]:
xs = np.linspace(-1000,1000,1000)
plt.plot(xs, max0(xs))
plt.plot(xs, LSE(xs))
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.legend(["$\max(0,x)$", "$LSE(x)$"])
plt.show()
Производная этой функции также важна и используется при обучении искуственных нейронных сетей.
$$\mathrm{softmax}(x) = LSE'(x) = \frac{e^x}{1+e^x}.$$График функции выглядит напоминает сглаженную ступеньку.
In [14]:
def softmax(x):
return np.exp(x)/(1+np.exp(x))
def step(x):
return 0.5*(np.sign(x)+1)
xs = np.linspace(-10,10,1000)
plt.plot(xs, step(xs))
plt.plot(xs, softmax(xs))
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.legend(["$step(x)$", "$softmax(x)$"])
plt.show()
На больших аргументах эту функцию снова не удается вычислить. Можно ли модифицировать softmax так, чтобы она работала на всем допустимом диапазоне для чисел с плавающей запятой?
In [15]:
xs = np.linspace(-1000,1000,1000)
plt.plot(xs, step(xs))
plt.plot(xs, softmax(xs))
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.legend(["$step(x)$", "$softmax(x)$"])
plt.show()
Попробуем оценить погрешность вычислений функции softmax. Пусть аргумент дан с машинной точностью, т.е. его отонсительная погрешность $\sigma x=\epsilon$.
Результат вычисления функции имеет вообще говоря погрешность отличную от погрешности аргумента. Например, для произвольной гладкой функции $f(x)$ и малой относительной погрешности $\sigma x$ аргумента $x$, предельная относительная погрешность результата оценивается по формуле:
$$\sigma f(x)\approx \kappa_f(x)\sigma x,$$где $\kappa_f(x)$ - число обусловленности для функции $f$ в точке $x$, вычисляемая по формуле:
$$\kappa_f(x)=\bigg|\frac{f'(x)x}{f(x)}\bigg|.$$Для экспоненты число обусловленности равно:
$$\kappa_{exp}(x)=\bigg|\frac{e^x x}{e^x}\bigg|=|x|.$$Значит погрешность вычисления числителя в $softmax$ оценивается так $$\sigma\,\mathrm{softmax}(x)\leq|x|\sigma x\leq|x|\epsilon.$$
Для оценки знаменятеля удобнее использовать абсолютную погрешность:
$$|\Delta(1+e^x)|\leq|\Delta 1|+|\Delta e^x|=0+e^x\sigma e^x\leq |x|e^x\epsilon,$$где мы воспользовались тем фактом, что погрешность константы равна нулю. Отсюда получаем предельную относительную погрешность:
$$\sigma(1+e^x)=\frac{|\Delta(1+e^x)|}{1+e^x}\leq \frac{e^x}{1+e^x}|x|\epsilon\leq |x|\epsilon.$$При делении предельные относительные погрешности складываются, следовательно:
$$\sigma\mathrm{softmax}(x)\leq \sigma e^x+\sigma(1+e^x)\leq 2|x|\epsilon.$$Результаты вычислений можно свести в таблицу:
| $a$ | $\Delta a$ | $\sigma a$ | ||||
|---|---|---|---|---|---|---|
| $x$ | $\epsilon$ | |||||
| $e^x$ | $e^x | x | \epsilon$ | $ | x | \epsilon$ |
| $1+e^x$ | $e^x | x | \epsilon$ | $ | x | \epsilon$ |
| $\frac{e^x}{1+e^x}$ | $2 | x | \epsilon$ |
Согласно нашим оценкам, погрешность вычислений softmax должна расти с увеличением аргумента. Попробуем оценить реальную погрешность. Функция softmax экспоненциально быстро стремиться к единице с ростом аргумента, т.е.
$$1-e^{-x} \leq \mathrm{softmax}(x) \leq 1,\quad x\to\infty.$$Т.е. для $x>40$ отличие softmax от 1 меньше $10^{-18}<\epsilon$, следовательно при оценки погрешности мы можем принять 1 за верное значение для достаточно больших $x$. Расчет на реальных данных показывает, что в пределах машинной точности отличий от 1 для больших аргументов нет.
In [16]:
xs = np.linspace(40, 700, 1000)
ys = softmax(xs)
err = np.abs((ys-1)/1) # относительная погрешность
print("Максимальная погрешность", np.max(err))
print("Минимальная погрешность", np.min(err))