Компьютерный практикум (магистратура)
Черновики лекций, весенний семестр 2022
http://mech.math.msu.su/~vvb/
Основная цель весеннего семестра 2022 г. --
система компьютерной математики SageMath
Подобные программные системы раньше назывались
системами компьютерной алгебры, поскольку изначально
подобные системы были ориентированы в большей степени на
символьные вычисления.
Наиболее известны Volfram Mathematika, MatLab, Maple,
Octave, Maxima, ...
Системы Volfram Mathematika, MatLab, Maple коммерческие,
также каждая из них использует собственный язык программирования
(на котором пишутся скрипты для решения конкретных задач).
Проект SageMath -- свободная система математики, которая
просто устанавливается под Linux. Особенности SageMath:
1) включает все существующие в мире свободные пакеты
(Singular для работы с многочленами, GAP для работы с группами
и т.д.), т.е. это оболочка для существующих свободных математических
пакетов;
2) в качестве языка программирования используется Python (сначала
это был Python2, c 2020 г. SageMath портирован на Python3);
3) SageMath даже не обязательно устаналивать на компьютере,
можно работать в облаке (даже с помощью смартфона).
Сайт: http://sagemath.org/
Мы начнем с языка Python3
Этот язык очень популярен в последнее время.
Почему?
Потому что он
1) очень легко учится;
2) имеет интерпретатор, можно сразу что-нибудь попробовать,
есть help, не обязательно читать книжки, есть tutorial
на сайте python.org. Python можно использовать как калькулятор;
3) python -- современный объектно-ориентированный язык программирования,
который включает некоторые возможности функциональных языков
программирования. (Функциональные языки -- это языки, в которых нет
переменных и циклов, главные объекты -- это функции и списки.);
4) python -- язык с динамической типизацией. Типы переменных и выражений
невозможно определить статически, по тексту программы (переменные
в Питоне вообще не имеют типов). В языке нет описаний типов;
5) большая часть приложений, которые используют машинное обучение,
нейросети и т.п., написаны на Python'е (на самом деле модули Питона
для машинного обучения написаны на C++, на Питоне пишут только самый
верхний уровень).
Мораль: если хочешь работать с искусственным интеллектом,
то Питона не избежать!
Как установить Python3
Следует просто установить интерпретатор Питона3
и затем можно установить какую-нибудь оболочку IDE,
для Питона (и для SageMath) лучше всего установить
Jupiter-lab
Первая программа на Питоне:
напечатать квадраты первых 20 натуральных чисел.
Списки (динамические массивы Питона)
smallPrimes = [2, 3, 5, 7, 11, 13, 17, 19]
Элементы списка доступна с помощью индекса:
idx = 0, 1, ..., len(smallPrimes) - 1
smallPrimes[4] = 11
Последний элемент списка:
smallPrimes[len(smallPrimes) - 1]
Если индекс отрицательный, то он отсчитывается от конца:
smallPrimes[-1] -- последний элемент списка
smallPrimes[-2] -- предпоследний элемент, ...
List comprehension (умное, или математическое задание списка)
Как задать в математике множество натуральных чисел
в отрезке 1..30, которые взаимно просты с 30 (т.е. не делятся
на 2, 3, 5)?
S = { x <- 1..30, gcd(x, 30) = 1 }
S = { 1, 7, 11, 13, 17, 19, 23, 29 }
В Питоне
s = [ x for x in range(1, 31) if math.gcd(x, 30) == 1 ]
ПЕРЕРЫВ до 10:50
Перерыв закончился.
List comprehension:
создаем список элементов.
Указываем выражение, зависящее от одного или нескольких
переменных. Для каждой переменной указываем
диапазон ее изменения +
можно указать фильтр, т.е. условие, при котором
объект добавляется в список.
Пример: список квадратов натуральных чисел от 1 до 20
s = [(x, x*x) for x in range(1, 21)]
Более сложный пример на использование List Comprehension
Получить список всех пифагоровых троек в диапазоне 1..100
Тройка (a, b, c), где a, b, c -- целые числа, пифагорова, если
1) a^2 + b^2 = c^2,
2) a < b,
3) a и b взаимно простые, т.е. gcd(a, b) = 1
s = [(a, b, c) for a in range(1, 101)
for b in range(a+1, 101)
for c in range(b+1, 101)
if math.gcd(a, b) == 1 and a*a + b*b == c*c]
Особенности списков в Питоне
1. Список -- это объект.
2. Списки можно изменять.
3. Отличие списков Питона от динамических массивов C++:
элементы списка могут иметь разные типы!
s = [1, "two", 3.0, [1, 1, 1, 1]]
Задача: получить список всех счастливых билетов длины n (четное число).
Номер билета счастливый, если сумма первых n/2 цифр билета равна
сумме последних n/2 цифр.
ВНИМАНИЕ!!!
В Python3 операция деления целых чисел выдает вещественный
результат.
>>> (2**20000)/(2**10000) == 2**10000
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
OverflowError: integer division result too large for a float
>>> (2**20000)//(2**10000) == 2**10000
True
Для целочисленного деления используется операция //
(она называется floordiv, отбрасывается дробная часть).
При работе с большими целыми числами надо обязательно
использовать операцию деления //, а не /
Задача: получить список всех счастливых билетов длины n (четное число).
Номер билета счастливый, если сумма первых n//2 цифр билета равна
сумме последних n//2 цифр.
Идея: у меня есть список целых чисел (цифр 0..9) длины n.
Я реализую функцию next(s), которая модифицирует этот список,
получая следующий список в лексикографическом порядке.
Фунция возвращает True, если это можно сделать, и False,
если этот список последний в лексикографическом порядке
[9, 9, ..., 9]
Первая задача:
написать функция, которая подсчитывает число счастливых билетов
длины n, n -- четное число, должна работать для
по крайней мере для всех n <= 16.
Домашнее задание
3 варианта, делаете тот вариант, остаток от деления на 3 от которого
совпадает с остатком от деления на 3 вашего номера по журналу.
1. Найти число счастливых билетов длины n <= 12.
def numLuckyTickets(n):
2. Проверить, является ли число полусовершенным.
def semiperfect(m):
Число совершенное, если оно равно сумме своих собственных
делителей.
6 == 1 + 2 + 3
24 == 1 + 2 + 3 + 6 + 12
Число полусовершенное, еслди оно равно сумме некоторого
подмножества множества своих собственных делителей.
Написать наивный алгоритм, основанный на решении задачи о рюкзаке.
Задача о рюкзаке:
дан список чисел и число n. Определить, можно ли из
некоторого подмножества этих чисел составить сумму, равную n.
3. Список s целых чисел содержит перестановку чисел 1, 2, ..., n.
Написать функцию nextPermutation(s), которая
преобразует тот же самый
список так, чтобы он содержал следующую в лексикографическом
порядке перестановку.
def nextPermutation(s):
. . .
Функция возвращает True, если это можно сделать, и False
в противном случае.
[1, 4, 5, 2, 6, 3] -->
[1, 2, 5, 3, 2, 6]
Для этой цели полезно:
s[i:j] -- подсписок, начиная с индекса i и заканчивая j-1
s[i:] -- от i до конца
s[:j] -- от начала до j-1
18 Feb 2022
Применение элементарной теории чисел в программировании
1. Повторим Python
Конструкция List Comprehension (появилась в Python'е
под влиянием функциональных языков программирования)
S = { x <- Z: 1 <= x <= 30, x не делится на 2, не делится на 3,
не делится на 5 }
s = [ x for x in range(1, 31) if x%2 != 0 and x%3 != 0 and x%5 != 0 ]
s -- множество обратимых элементов кольца Z_30
Эта конструкция использует более общую конструкцию
generator expression
Конструкция generator expression на самом деле использует понятие
generator function
Конструкция generator expression сводится к объекту типа iterator
У такого объекта долны быть реализованы два метода __iter__()
и __next__(), которые вызываются с помощью функций iter(i) и
next(i).
class Abcd:
. . .
def __iter__(self):
return self
def __next__(self):
. . .
# Возвращает очередной объект
# По окончанию возбуждает исключение типа
# StopIteration
Генераторная функция
Это функция, которая выдает некоторую последовательность
элементов, по одному элементу за один раз.
Функция вместо того, чтобы вернуть результат с помощь return,
выдает очередной объект с помощью yield
Эта функция обычно используется в цикле for, т.е. по ней
стороится объект типа iterator; этот объект запоминает все локальные
переменные (Frame) функции и точку, на которой оа прервалась
при выполнении yield; при последующем обращении функция возобновляет
работу с прерванного места.
Пример:
напишем генераторную функцию, которая перебирает целочисленные
точки первого квадранта плоскости
def intPoints(n):
(0, 0),
(1, 0), (0, 1),
(2, 0), (1, 1), (0, 2),
(3, 0), (2, 1), (1, 2), (0, 3),
. . .
... (0, n)
Эту функцию мы используем для нахождения хорошего приближения
числа pi в виде m/n, |pi - m/n| <= 1/n^2
--------------------------------------------------
import math
def intPoint(n = 1000):
for k in range(n + 1):
x = k
y = 0
yield (x, y)
while x > 0:
x -= 1
y += 1
yield (x, y)
def piApproximation(eps = 0.00001):
'''Compute a fraction m/n so that
abs(pi - m/n) <= eps and
abs(pi - m/n) <= 1/(n*n)'''
for (m, n) in intPoint(1000000):
if (n == 0):
continue
err = abs(math.pi - m/n)
if err <= eps and err <= 1/(n*n):
return (m, n)
return None
--------------------------------------------------
Вот результат работы этой программы:
>>> from piApprox import *
>>> piApproximation()
(355, 113)
>>> 355/113
3.1415929203539825
>>> math.pi
3.141592653589793
--------------------------------------------------
Теория чисел
1. Основной объект -- кольцо вычетов по модулю m
Z/mZ = Z_m
Это фактор множество кольца целых чисел по отношению
эквивалентности
x == y <===> m | (x - y)
def
Элементами фактор-множества является классы эквивалентности
В программировании мы как бы работаем с целыми числами,
учитывая, что числа x и y представляют один и тот же
элемент Z_m, когда m | (x - y)
Всегда можно x --> x%m
Если в каждом классе эквивалентности выбрать один элемент,
то получим систему остатков. Систем остатков бесконечно много,
но обычно рассматривают только две системы:
1) неотрицательная
Пример: Z_5 { 0, 1, 2, 3, 4 }
2) симметричная система
{ -2, -1, 0, 1, 2 }
Если m четное, она не совсем симметрична
Например, при m = 10 это
{ -5, 4, 3, 2, 1, 0, 1, 2, 3, 4 }
Традиционно отрицательных остатков на 1 больше, чем
положительных
d = m/2
d + d == 0 (mod m) d = -d
В частности, для чисел типа int в языке C
2^31 -- это максимальное по модулю
отрицательное число -2147483648
2^31 - 1 -- это максимальное по модулю
положительное число 2147483647.
int в С == Z_m где m = 2^32 = 4294967296.
R -- ассоциативное кольцо I -- идеал в кольце
R/I -- фактор-кольцо
x ~ y <==> (x - y) <- I
Кольцо Z является кольцом главных идеалов:
Любой идеал в кольце Z равен m*Z
В Z все фактор-кольца имеют вид m*Z
Z_m = Z/m*Z
В программировании практически никогда не работают с кольцом
целых чисел. Всегда работают в кольца Z_m.
Например, "целые" числа типа int в языке C -- это элементы
фактор-кольца Z_m при m = 2^32
Предупреждение:
как возвести k в степень n по модулю m?
Ни в коем случае нельзя писать так, как пишет половина
студентов:
p = k**n % m
В криптографии n -- 400-значное целое число. Число k**n
невозможно записать ни в какой памяти!
В Питоне есть встроенная функция pow с тремя аргументами:
p = pow(k, n, m)
В программировании есть несколько базовых алгоритмов,
на которых основана любая работа в кольцах вычетов
1. Быстрое возведение в степень по модулю m (в кольце Z_m)
Дано: x, n >= 0, m <-Z
Надо: вычислить x^n в кольце Z_m
Применим схему посроения цикла с помощью инварианта:
p * y^k == x^n (mod m)
Начальные значения:
p = 1, y = x, k = n (инвариант выполняется)
Условие завершения цикла:
k == 0, тогда ответ в p
Действие в теле цикла, уменьшающее k и при этом сохраняющее
инвариант:
k четное => k /= 2, y *= y
k нечетное => k -= 1, p *= y
Расширенный алгоритм Евклида
============================
Теорема:
Пусть m, n <- Z, m != 0 или n != 0 (пара (0, 0) исключается).
Тогда d = gcd(m, n) выражается в виде линейной комбинации
чисел m и n с целыми клоеффициентами:
d = u*m + v*n,
где u, v <- Z
Расширенный алгоритм Евклида позволяет по числам m, n
вычислить тройку чисел (d, u, v)
ПЕРЕРЫВ 10 минут до 10:40
Перерыв закончился.
Реализуем Расширенный алг. Евклида, применив схему построения
цикла с помощью инварианта.
Идея: применяем обычный алгоритм Евклида
(a, b) --> (b, r), где r -- остаток от деления a на b
В начале (a, b) = (m, n)
В любой момент времени храним коэффициенты выражения a и b
через m и n.
Инвариант цикла:
gcd(a, b) = gcd(m, n),
a = u1*m + v1*n,
b = u2*m + v2*n
Начальные значения:
(a, b) = (m, n)
u1 = 1, v1 = 0,
u2 = 0, v2 = 1
Действия в цикле, сохраняющие инвариант
(a, b) --> (b, r), где a = q*b + r
Тогда r = a - q*b = u1*m + v1*n - q*(u2*m + v2*n) =
= (u1 - q*u2)*m + (v1 - q*v2)*n
(u1, u2) --> (u2, u1 - q*u2)
(v1, v2) --> (v2, v1 - q*v2)
Значение расширенного алгоритма Евклида в том,
что он позволяет вычислить
обратный элемент к n в кольце вычетов Z_m
Элемент n обратим в Z_m <==> когда gcd(n, m) = 1
1 = u*m + v*n Так как u*m == 0 (mod m), то
1 == v*n (mod m), то есть v -- обратный элемент к n в Z_m
Китайская теорема об остатках
Пусть m = m1*m2*...*m_k, где m_i, m_j взаимно простые при i != j
gcd(m_i, m_j) = 1 при i != j
Тогда кольцо Z_m изоморфно прямой сумме колец Z_mi, i = 1, ..., k
~
Z_m = Z_m1 + Z_m2 + ... + Z_mk
То есть элемент Z_m, т.е. число x, представляется в виде строки
x == (x1, x2, ..., xk)
x_i == x (mod m_i)
Операции в прямой сумме выполняются покомпонентно.
Доказательство:
Рассмотрим гомоморфизм
f: Z --> Z_m1 + Z_m2 + ... + Z_mk
x --> (x1, x2, ..., xk)
Ядро этого гомоморфизма -- это элементы x,
которые делятся на m1, m2, ..., mk
Поскольку они взаимно простые, то x делится на m1*m2*...*mk = m
Значит, ker f = m*Z
По теореме о гомоморфизме
Z/Ker(f) == Im(f) -- инъекция.
|Z/Ker(f)| = m и |Z_m1 + Z_m2 + ... + Z_mk| = m
=> f -- сюрьективно => f биективно
Доказательство не конструктивное, но есть алгоритм, который
для строки (x1, x2, ..., x_k) вычисляет x:
x == x_i (mod m_i)
Смысл китайской теоремы об остатках состоит в том,
что мы сводим изучение кольца Z_m для произвольного m
к изучению примарных колец вида
Z_p^l
где p -- простое число.
Любой число m представимо в виде произведения степеней простых
m = p1^e1 * p2^e2 * ... * pk^ek
Z_m = Z_p1^e1 + ... + Z_pk^eK
Пример: сколько корней из 1 в кольце Z_105 ?
Подсказка: 105 = 3*5*7
x*x == 1 (mod 105)
Z_105 = Z_3 + Z_5 + Z_7
x = (x1, x2, x3)
x^2 = (x1^2, x2^2, x3^2)
x^2 == 1 (mod 105) <==> x_i^2 == 1
Сколько квадратных корней из 1 в любом поле?
x^2 - 1 = 0
Ответ: 1, -1
Значит x^2 == 1 => x = (+-1, +-1, +-1)
Ответ: 8 элементов
(1, 1, 1)
(-1, 1, 1)
(1, -1, 1)
. . .
(-1, -1, -1)
Алгоритм для Китайской теоремы об остатках
Ограничимся случаем k = 2
Дано: числа r1, r2 -- любые
числа m1, m2 -- взаимно простые gcd(m1, m2) = 1
Надо: найти x:
x == r1 (mod m1)
x == r2 (mod m2)
Решение:
Из первого сравнения вытекает, что
x = r1 + s*m1, s -- любой целое число
Подберем s так, чтобы выполнялось второе сравнение
r1 + s*m1 == r2 (mod m2)
s*m1 == r2 - r1 (mod m2)
Элемент m1 обратим в кольце Z_m2
Найдем v: v*m1 == 1 (mod m2)
Домножим равенство на v
v*s*m1 == v*(r2 - r1) (mod m2)
s == v*(r2 - r1) (mod m2)
Напишем это на Питоне:
def chineseReminderAlg(r1, r2, m1, m2):
v = invmod(m1, m2)
s = (v*(r2 - r1))%m2
return r1 + s*m1
Замечательные теоремы элементарной теории чисел
Малая теорема Ферма
Пусть m -- простое число. Тогда для всякого b != 0 (mod m)
b^(m-1) == 1 (mod m)
Частный случай теоремы Лагранжа:
для всякой конечной группы G порядка n
для всякого элемента x <- G
x^n = 1
Порядок любого элемента делит порядок группы.
В Малой теореме Ферма рассматривается группа ненулевых элементов
кольца Z_m по умножению при простом m. Ее порядок равен m-1.
Другие формулировки Малой теоремы Ферма
Для всякого b и для простого p
b^p == b (mod p)
или в чуть более общем виде
b^(k*(p-1) + 1) == b (mod p) для любого b
Именно эта формулировка Малой теоремы Ферма используется
в схеме RSA кодирования с открытым ключом.
Обобщение Малой теоремы Ферма
Функция Эйлера phi(m) = | U_m |
где U_m -- группа обратимых элементов кольца Z_m
Выведем формулу для функции Эйлера
m = p1^e1 * p2^e2 * ... * p_k^e_k
Тогда
phi(m) = (p1 - 1)*p1^(e1-1) * ... * (p_k - 1)*p_k^(e_k - 1)
В частном случае
m = p1 * p2 * p3 * ... * p_k
m свободно от квадратов. Тогда
phi(m) = (p1 - 1)(p2 - 1)...(p_k - 1)
Теорема Эйлера
Пусть b взаимно просто с m. Тогда
b^phi(m) == 1 (mod m)
Частный случай теоремы Эйлера
Пусть m свободно от квадратов,
m = p1 * p2 * ... * p_k, p_i простые.
Тогда для любого b <- Z и любого s <- Z, s >= 0
b^(s*phi(m) + 1) == b (mod m)
Возведение в степень s*phi(m) + 1 является
тождественным отображением в кольце Z_m.
Доказательство:
следствие Малой теоремы Ферма в третьей формулировке
и Китайской теоремы об остатках.
Вывод формулы Эйлера
По Китайской теореме об остатках, достаточно доказать,
что phi(p^e) = (p - 1)*p^(e - 1)
(достаточно доказать ф-лу Эйлера для примарного кольца).
Сколько обратимых элементов в Z_p^e ?
Проще сосчитать, сколько необратимых элементов.
Необратимы элементы, которые делятся на p,
т.е. каждый p-й элемет.
Всего их (p^e/p) = p^(e-1)
Число обратимых = число всех - число необратимых =
p^e - p^(e-1) = (p - 1)p^(e - 1)
Схема кодирования с открытым ключом (асимметричное кодирование)
Классические шифры
1. Шифр Цезаря -- циклический сдвиг на 3 позиции по алфавиту
ABCDEF...Z
DEFGHI...C
Для расшифровки надо было сделать обратный сдвиг
2. Шифр Цезаря -- частный случай подстановочного шифра
Алфавит S
Шифр -- биекция f: S --> S
Расшифровка f^(-1): S --> S
Любой подстановочный шифр легко взламывается с помощью частотного
анализа
3. Шифр Виженера -- с перекрытием текста.
Взламывается с помощью методов Касицкого (основанного на
биграммах) и Фридмана (основан на индексе совпадения).
4. Роторные шифровальные машины (они применялись уже во время
первой мировой войны). Фашистская Германия использовала
шифровальную машину Enigma.
Это была электоро-механическая машина -- внешне пишушая машинка,
шифровали двое. Первый человек нажимал букву на клавиатуре,
загоралась лампочка, соответствующая зашифрованной букве,
второй человек ее записывал. Расшифровка была такой же.
Устройство: 3 ротора + рефлектор
Каждый ротор осуществлял отображение (перестановку)
конечного множества из 26 элементов --
26 контактов слева и справа, внутри соединенных
проводами.
s1 -- перестановка 1-го ротора
s2 -- перестановка 2-го ротора
s3 -- перестановка 3-го ротора
t -- перестановка рефлектора
t -- произведение независимых транспозиций
t^2 = 1 идемпотент
Первый ротор поворачивался на 1 позицию после каждой буквы
Второй ротор поворачивался после полного поворота
первого ротора (через 26 букв)
Третий ротор поворачивался после полного поворота
второго ротора (через 26*26 букв)
Ключ шифрования -- начальное положение роторов.
На каждом шаге выполнялась перестановка
sigma = s1 s2 s3 t s3^-1 s2^-1 s1^-1
sigma^2 = 1
Для расшифровки надо было снова пропустить зашифрованное
сообщение через машину, восстановив начальное положение роторов.
Начальное положение роторов было указано в шифровальной книге
для каждого дня года.
Enigma была взломана польскими математиками под руководством
Мариана Решевского. После оккупации Польши работы были переданы
Англии.
============================================================
25 Feb 2022
Применение теории чисел в криптографии
История Энигмы и ее расшифровки убеждает в несостоятельности
классических шифров, основанных на заменах букв алфавита.
Первая идея в современных шифрах:
блочное шифрование. Сообщение сначала представляется в виде
последовательности нулей и единиц (битов), затем эта последовательность
разбивается на блоки фиксированного размера (например, 128 битов).
Каждый блок шифруется отдельно (либо независимо от предыдущего блока,
либо есть схемы, где образ блока вычисляется с использованием также и
предыдущего блока). Получается зашифрованный блок того же размера.
Такого рода шифр, например, DES (Data Encryption Standard),
использует набор преобразований битов, 16 раундов. Каждое преобразование
обратимо, преобразования задаются ключом. Идея этих преобразований
определяется двумя словами: запутывание и перемешиание (confusing and
diffusing).
Симметричные шифры были придуманы в 50-60 гг. и обычно имели статус
государственного секрета. Но в начале 70-х годов в MIT была предложена
схема асимметричного кодирования, основанная на простых результатах
теории чисел: RSA (Rumley, Shamir, Adleman).
Она называется "Кодирование с открытым ключом" (Public Key Encryption).
Она никогда не имела статуса государственного секрета.
В схеме RSA используются два ключа. Открытый ключ (public key) известен всем
и позволяет кому угодно зашифровать и послать сообщение.
Но расшифровать его может только обладатель секретного ключа (private,
or secret key).
Главная идея RSA -- блок битов интерпретировать как двоичную запись
целого числа, вернее, кольца Z_m вычетов по модулю m
(фактор кольцо кольца целых чисел Z по главному идеалу mZ).
Кольцо Z_m конечно и в нем можно быстро выполнять любые операции,
в часности, возведение в степень.
Асимметричное кодирование:
Кодирующая процедура -- это инъективное отображение
E: Z_m --> Z_k, k >= m (в RSA E: Z_m --> Z_m -- биективное отображение)
Отображение E задано открытым ключом и известно всем.
Обратное отображение
D: Z_k --> Z_m
является секретным и легко не восстанавливается по E
(секретный ключ современными алгоритмами невозможно вычислить по
открытому ключу).
D E = 1 (тождественное отображение Z_m --> Z_m)
Для всякого x
y = E(x) -- зашифрованное сообщение
D(y) = D(E(x)) = x
В RSA
E: Z_m --> Z_m
D: Z_m --> Z_m
Не только DE = 1, но и ED = 1
Это позволяет применять шифрование в обратном направлении,
в частности, удостоверять отправителя сообщения.
x -- открытое сообщение
Обладатель секретного ключа его шифрует секретнолй прцедурой D:
y = D(x)
Любой может его прочесть, поскольку процедура E открыта
x = E(y)
Но послать его может только обладатель секретной процедуры D.
На той же идее основана электронная подпись:
Посылаем открыто файл T и к нему электронную подпись:
вычисляются
1) хэш-функция от содержимого T,
2) длина T,
3) номер транзакции,
4) текущие дата и время и т.п.
и к этим данным применяется секретная процедура D.
Любой может восстановить эти данные, применив открытую процедуру E
и проверить, сопадают ли они с данными подписанного файла T.
Если нет, то файл подменен мошенниками.
Такие же идеи используются в криптовалютах и т.п.
Перейдем к описанию схемы шифрования RSA.
Рассмотрим 2 больших простых числа p, q (порядка 200 десятичных цифр каждое).
Пусть m = p*q
m открыто, разложение m на множители секретно
Взлом ключа сводится к разложению m на множители, которое
невозможно сделать современными алгоритмами
(пока не реализован квантовый компьютер!).
Безем любое число e <- Z_phi(m), обратимое в кольце Z_phi(m)
Здесь phi(m) -- функция Эйлера
phi(m) = (p - 1)*(q - 1)
Вычислим с помощью расширенного алгоритма Евклида обратный элемент d к e
по модулю phi(m):
d*e == 1 (mod phi(m))
Шифрующее отображение:
E: x --> x^e (mod m)
Дешифровка:
D: y --> y^d (mod m)
Покажем, что это взаимно обратные отображения
E(x) = y == x^e (mod m)
D(y) == y^d == (x^e)^d == x^(e*d) (mod m)
Имеем
d*e == 1 (mod phi(m))
d*e = 1 + k*phi(m)
По теореме Эйлера (вариант для числа m, свободного от квадратов)
для всякого x
x^(1 + k*phi(m)) == x (mod m)
Мы доказали, что для всякого x
D(E(x)) = x
Значит, DE = 1 и ED = 1
Предложение доказано.
Итак, открытый ключ -- это пара (m, e)
E: x --> x^e (mod m)
Секретный ключ -- это пара (m, d)
D: y --> y^d (mod m)
Пример:
>>> p = 986064567478133119510880837311
>>> q = 1241326039840269783519222613267
>>> m
12240276245744394650079821223327782866261146467026504972050
>>> m = p*q
>>> phi = (p - 1)*(q - 1)
>>> phi
1224027624574439465007982122330550896018796243799620393754460
>>> e = 1153932867384633370819307154242378038591183291
>>> d = invmod(e, phi)
>>> d
375184827923908186491134676097352998468766988827636275736371
>>> d*e % phi
1
Открытый ключ -- пара
m = 12240276245744394650079821223327782866261146467026504972050
e = 1153932867384633370819307154242378038591183291
Секретный ключ -- пара
m = 12240276245744394650079821223327782866261146467026504972050
d = 375184827923908186491134676097352998468766988827636275736371
Проверим шифрование и расшифровку сообщения
>>> x = 12345678910111213141516
Шифруем:
>>> y = powmod(x, e, m)
>>> y
790416752955372286760677863659087748274935561825490735814409
Это зашифрованное сообщение.
Расшифровка:
>>> z = powmod(y, d, m)
>>> z
12345678910111213141516
>>> z == x
True
Получили исходное сообщение.
ПЕРЕРЫВ 10 минут до 10:50
Какие задачи возникают в связи со схемой RSA?
Чтобы сгенерировать пару ключей, надо уметь генерировать
большие простые числа p, q (причем не любые простые числа подходят,
надо, чтобы m = p*q невозможно было разложить на множители
существующими алгоритмами, в частностности, p-1 и q-1 не должны
быть гладкими (smooth) числами).
Простых чисел очень много: формула Чебышева
(число простых <= n) ~ n/ln(n)
Грубо говоря, в отрезке длины log(n) в среднем содержится одно
простое число.
Значит, если у нас есть тест простоты (функция isPrime(m) == True/False),
то мы можем быстро найти случайное простое число заданной длины.
Первая задача: реализовать тест простоты.
Вторая задача: для взлома ключа надо вычислить функцию Эйлера
phi(m) = (p - 1)*(q - 1)
(чтобы вычислить обратный к e элемент d в кольце Z_phi(m)),
для этого надо разложить m на множители:
m = p*q
Обратно, зная функцию Эйлера, легко можно разложить число m на
множители (есть много разных способов).
Таким образом, задача факторизации эквивалентна задаче вычисления
функции Эйлера.
Итак, для взлома ключа надо решать задачу факторизации целого числа
(факторизация == разложение на множители).
Перейдем к тестам простоты.
Первый "тест" -- Малая теорема Ферма
m -- простое, то для всякого b != 0 (mod m)
b^(m-1) == 1 (mod m)
Это необходимое условие простоты.
Довольно долго существовало заблуждение,
что это и достаточное условие, причем даже только
для основания 2 (так считали древние греки):
Если 2^(m-1) == 1 (mod m) ==> m простое
Это неверно! Впервые контрпример привел французский математик П.Сарруз
в середине XIX века:
m = 341 = 31*11 -- не простое число
2^(m-1) (mod m) == 2^340 = (2^10)^34 = 1024^34 (mod 341) ==
1024 = 341*3 + 1 == 1 (mod 341)
== (341 + 1)^34 (mod 341) == 1^34 == 1 (mod 341)
Число 341 веден себя как простое число в Малой теореме Ферма,
если основание степени равно 2.
Но для основания 3 Малая теорема Ферма уже не выполняется:
>>> powmod(3, 340, 341)
56
Но в 1910 американский математик Кармайкл нашел удивительное число
561, которое не яляется простым, но ведет себя как простое в Малой
теореме Ферма. В его честь такие числа называются Кармайкловыми.
Число m называется кармайкловым (carmichael numbers), если
для всякого b, взаимно простого с m, gcd(b, m) = 1,
выполняется утверждение Малой теоремы Ферма:
b^(m-1) == 1 (mod m)
Мораль: для кармайкловых чисел Малая теорема Ферма не годится
в качестве теста простоты.
Кармайкловы числа:
561, 1105, 1729, ...
их бесконечно много.
Они устроены очень просто:
Теорема:
нечетное число m кармайклово <==> когда выполняются все 3 условия:
1) m свободно от квадратов
m = p1*p2*...*p_k,
2) k >= 3,
3) для всякого i выполняется
(p_i - 1) | (m - 1)
Доказательство достаточности сразу получается из Китайской теоремы
об остатках
b^(m-1) == 1 (mod m)
b == (b1, b2, ..., b_k)
b^(m-1) == (b1^(m-1), b2^(m-1), ..., b_k^(m-1))
b_i^(m-1) == 1, поскольку (p_i - 1) | (m - 1)
m - 1 = (p_i - 1)*s
b_i^(m-1) = (b_i^(p_i - 1))^s == (1)^s == 1 (mod p_i)
Пример: 561 = 3*11*17
(3-1) | 560
(11-1) | 560
(17-1) | 560, поскольку 560 = 35*16
1729 = 7*13*19
6 | 1728,
12 | 1728,
18 | 1728
Но, оказывается, тест Ферма можно чуть-чуть подправить так,
что получится вероятностный тест простоты Миллера--Рабина.
Что такое вероятностный тест?
В тесте Миллера-Рабина проверятся простота числа m.
Тест использует датчик случайных чисел, с помощью которого мы
генерируем целое число b, которое используется как основание
степени в кольце Z_m.
Тест для фиксированного b выдает один из двух ответов:
1) число m составное, при этом приводится доказательство этого факта;
2) не знаю.
Причем для составного числа m вероятность второго ответа не превышает 1/4.
Главное, что эта оценка сверху вероятности не зависит от m.
Пусть m -- составное. Пусть мы применили 20 раз тест и получили 20
ответов "не знаю". Тогда вероятность такого события
не превосходит (1/4)^20 = 2^(-40), почти что ноль. Делая большее
число тестов, эту вероятность можно сделать меньше любого eps.
Значит, скорее всего в этом случае число простое, но, увы,
мы не получаем доказательства этого факта.
Но, тем не менее, такие "псевдопростые" числа достаточны для криптографии.
Отметим, что в 2006 г. был найден детерминированный метод доказательства
простоты числа m, который работает за полиномиальное время в зависимости
от числа цифр в записи m. (Оценка порядка O(n^6).)
Алгоритм Миллера--Рабина.
Проверяем на простоту нечетное число m.
Берем случайное число b <- Z_m
Малая теорема Ферма: b^(m-1) == 1 (mod m)
Вычисляем число b^(m-1) следующим образом:
поскольку m-1 -- четное число, то
m - 1 = 2^s * t, где t -- нечетное число.
Вычиcляем следующую последовательность в кольце Z_m:
x0 == b^t, x1 == x0^2, x2 == x1^2, ... x_s = x_s-1^2 (mod m)
Имеем
x_s == b^(m-1) (mod m)
Рассмотрим последовательность, каждый следующий элемент
в ней есть квадрат предыдущего (mod m):
x0, x1, x2, ..., x_s
Утверждение:
если эта последовательность кончается не единицей
*, *, *, ..., *
или если она содержит фрагмент вида
*, *, *, ..., *, 1, 1, ..., 1
где * != +-1 (mod m)
то число составное, тест выдает ответ "составное".
Ответ "не знаю" выдается в случаях, когда последовательность имеет
один из двух видов:
1, 1, 1, ..., 1
или
*, *, *, ..., *, -1, 1, ..., 1
В этом случае остается подозрение, что число m простое
(вероятность подобной последовательности для составного m
не превышает 1/4).
Покажем, что тест никогда не врет, т.е. если выдается ответ "составное",
то число m на самом деле составное.
Действительно, если последовательность кончается не единицей
*, *, *, ..., *
то не выполняется Малая теорема Ферма;
если же последовательность имеет вид
*, *, *, ..., *, 1, 1, ..., 1
x0, ... x_i, 1, ..., 1
где x_i != +- 1, но x_i^2 == 1 (mod m),
то число x_i является квадратным корнем из 1.
Но в случае простого m кольцо Z_m является полем, а в любом поле нет
других квадратных корней из 1, кроме 1 и -1.
Почему? Корни яляются решениями уравнения
x^2 - 1 = 0,
а по теореме Безу число корней полинома не превосходит его степени.
Вторая часть утверждения -- что для составного m вероятность
последовательности
1, 1, 1, ..., 1
или
*, *, *, ..., *, -1, 1, ..., 1
не превышает 1/2, вытекает из того, что число элементов
всякой СОБСТВЕННОЙ подгруппы группы U не превосходит |U|/2.
Таким образом доказывается неравенство с константой 1/2
(подробности доказательства мы не приводим).
Немного более сложные рассуждения позволяют получить константу 1/4
(что не принципиально).
Следующий раз рассмотрим алгоритмы факторизации
целых чисел и алгоритм извлечения квадратного корня в поле
вычетов по модулю p, где p -- большое простое число.
Ревков
Просьба в теме письма писать слово Магистратура.
====================================================
4 Mar 2022
Применение теории чисел в программировании
Была рассмотрена схема кодирования с открытым ключом RSA
(математики из MIT).
m = p*q, где p, q -- дольшие простые числа
m открыто, p, q -- секретные
e -- любой обратимый элемент в кольце Z_phi(m)
phi(m) = (p - 1)*(q - 1) -- функция Эйлера
d -- обратный к e элемент
d*e == 1 (mod phi(m))
e открыто, d секретно
Отображение кодирования (Encoding)
E: x --> x^e (mod m)
Обратное отображения (Decoding)
D: y --> y^d (mod m)
Эти отображения взаимно обратные:
D*E = E*D = 1 (тождественное отображение Z_m --> Z_m)
x -- исходное сообщение
y == x^e (mod m) -- зашифрованное сообщение
z == y^d (mod m) == x -- расшифровка
Открытый ключ (public key): пара (m, e)
Секретный ключ (private key): пара (m, d)
Для генерации ключа надо сгенерировать два больших простых числа p, q
Затем вычисляем m = p*q
phi(m) = (p - 1)*(q - 1)
Берем случайное e, смотрим, обратимо оно или нет: gcd(phi(m), e) = 1
С помощью расширенного алгоритма Евклида вычисляем d = inverse(e) (mod phi(m))
Для взлома открытого ключа (m, e) надо вычислить d, для этого надо
вычислить ф-ю Эйлера, а для этого в свою очередь надо разложить m
на множители -- а это сложная задача.
Насколько сложная, пока никто не знает (например, мы не знаем,
является ли эта задача NP-полной).
_ _
(x*y*z*t*s) + (.......) + .... + ( .... ) Дизъюнктивная нормальная форма
_ _
(x+y+z+t+s) * (.......) * .... * ( .... ) Конъюнктивная нормальная форма
Конкурс SAT-решателей
(a SAT solver is a computer program which aims to solve
the Boolean satisfiability problem. On input a formula over
Boolean variables, such as "(x or y) and (x or not y)",
a SAT solver outputs whether the formula is satisfiable,
meaning that there are possible values of x and y
which make the formula true, or unsatisfiable, meaning
that there are no such values of x and y.)
С помощью SAT-решателя была
решена отрицательно задача о хроматическом числе плоскости
для для числа цветов n = 4: существует плоский граф, длины
ребер которого равны единице, вершины которого нельзя раскрасить
в 4 цвета так, чтобы смежные вершины имели разные цвета.
Для n=5 и n=6 проблема остается пока открытой.
На данный момент есть алгоритмы, которые факторизуют 100-значные числа
(может, 150-значные) с субэкспоненциальным временем работы.
Уж точно 400-значное десятичное число невозможно разложить на множители,
но сгенерировать 200-значные простые числа несложно.
Мораль: надо изучить алгоритмы факторизации.
Это задача интересна и помимо практической криптографии, просто для
математики.
Перерыв 10 минут до 11:05
Алгоритмы факторизации
Совсем простой -- колесо со спицами (совсем слабый)
Простые, на классные два алгоритма Полларда (методы Монте-Карло:
если повезет, то разложим число, но это не гарантировано):
1) rho-алгоритм Полларда;
2) p-1-алгоритм Полларда.
Усиление p-1-алгоритма Полларда с помощью использования
эллиптических кривых над конечным полем или кольцом:
алгоритм факторизации Ленстра.
Алгорит "квадратичное решето" (он был лучшим в конце 80-х годов).
Сейчас первенство удерживает алгоритм "квадратичное решето над
конечным расширением (квадратичным?) поля Z_p.
1. rho-алгоритм факторизации целого числа Полларда,
основан на поиске цикла в рекуррентной последовательности.
Дано: число m > 0, m <- Z
Надо: найти нетривиальный делитель числа m
d | m, 1 < d < m
Рассмотрим полиномиальное отображение
f: Z_m --> Z_m
f: x |--> x^2 + 1 (mod m)
Берем случайный элемент x0 <- Z_m и строим последовательность
x0, x1 = f(x0), x2 = f(x1), x3 = f(x2), ..., x_i = f(x_i-1), ...
x0, x1, x2, x3, x4, ...
Поскольку |Z_m| < бесконечности, то последовательность {x_i} циклическая.
Найдется n, d: x_n = x_n+d
d -- длина цикла.
(Следовательно, для всех s >= n x_s == x_s+d (mod m)
Пусть p | m (мы p пока не знаем). Рассмотрим ту же последовательность,
но по модулю p (это корректно, т.к. p | m):
_ _ _ _ _
x0, x1, x2, x3, x4, ...
_
x_i == x_i (mod p) (это образы элементов x_i при каноническом
эпиморфизме Z_m --> Z_p).
_
Последовательность {x_i} тоже циклическая, причем, скорее всего,
ее цикл короче: найдутся k, t: x_k == x_k+1 (mod p)
Имеем x_i == x_i+t (mod p)
x_i != x_i+t (mod m)
Обозначим x_i = y, x_i+t = z
p | (y - z)
m не делит (y - z) ==>
p | gcd(m, y-z)
m не делит gcd(m, y-z) ==> d = gcd(m, y-z) -- это нетривиальный
делитель числа m
Мы свели задачу к поиску цикла в рекуррентной последовательности.
Самая простая схема:
сравниваем
x0 <-->x1
x1 <-->x3
x2 <-->x5
x3 <-->x7
x4 <-->x8
. . .
Сравнение y и z -- это вычисление gcd(m, y-z)
Обычный примитивный алгоритм факторизации или проверки простоты:
Дано: m
Надо: найти нетривиальный делитель d | m
Проверяем, делится ли m на 2
Если нет, то делим на 3, 5, 7, 9, 11, 13, ... до квадратного корня
из m.
Пусть n = 2*3*5 = 30
Выпишем все числа от 1 до 29, которые не делятся на 2, 3, 5,
т.е. взаимно просты с 30:
dividers = [ d for d in range(1, 30) if gcd(d, 30) == 1 ]
[1, 7, 11, 13, 17, 19, 23, 29]
Сначала делим на 2, 3, 5
Если m не делится на 2, 3, 5, то выполняем пробные деления на числа
7, 11, 13, 17, 19, 23, 29,
31, 37, 41, 43, 47, 49, 53, 59,
61, 67, 71, 73, 77, 79, 83, 89,
. . .
k*30+1, k*30+7, k*30+11, ..., k*30+29,
. . .
Всего доля пробных делителей 8/30 (меньше, чем 1/2 в наивном алгоритме)
===============================================
11 Mar 2022
Была рассмотрена схема кодирования с открытым ключом RSA
m = p*q
phi(m) = (p-1)(q-1)
e -- обратимо в Z_phi(m)
d -- обратный к e элемент
e*d == 1 (mod phi(m))
Encyption E: Z_m --> Z_m
E: x --> x^e (mod m)
Decription D: Z_m --> Z_m
D: y --> y^d (mod m)
DE = ED = 1 -- тождественное отображение в Z_m
Был рассказан вероятностный тест простоты числа, с помощью которого
можно находить большие простые числа => мы легко можем генерировать
пары ключей.
Задача взлома ключа сводится к вычислению функции Эйлера phi(m),
что эквивалентно задаче факторизации числа m.
Знаем разложение на множители числа m
m = p*q ==>
знаем функцию Эйлера
phi(m) = (p - 1)(q - 1)
Обратно, если знаем функцию Эйлера, то можно разложить m на множители:
(p - 1)(q - 1) = phi
p*q = m
Система из двух уравнений с двумя неизвестными, которая легко решается
(сводится к квадратному уравнению).
Задача факторизации
===================
Дано число m, про которое известно, что оно составное.
Найти нетривиальный делитель d | m
Есть несколько красивых алгоритмов факторизации.
Прошлый раз мы рассмотрели
1) rho-алгоритм Полларда (метод "Монте-Карло"), основанный
на поиске цикла в последовательности
f: Z_m --> Z_m
f: x --> x^2 + 1 (mod m)
x0, x1, x2, x3, ... <- Z_m
x_i+1 = f(x_i)
Для p | m, то можно рассмотреть ту же последовательность, но по (mod p)
_ _ _ _
x0, x1, x2, x3, ... <- Z_m
_
x_i == x_i (mod p) Z_m --> Z_p
_
Последовательность { x_i }, скорее всего, имеет цикл длины, меньшей,
чем последовательность { x_i }
x_i == x_i+l (mod p)
x_i != x_i+l (mod m)
p | (x_i - x_i+l)
m не делит (x_i - x_i+l)
=> p | gcd(x_i - x_i+l, m)
m не делит gcd(x_i - x_i+l, m)
=> d = gcd(x_i - x_i+l, m) -- нетривиальный делитель m.
Оценка скорости работы.
Она основана на оценке длины цикла в орбите произвольного
элемента x0 <- Z_m при действии полиномиального отображения
f: x --> x*x + 1 (mod m)
Известно, что для полиномиального отображения Z_m -- > Z_m
длина цикла в орбите произвольного элемента имеет порядок O(m^(1/2)).
Наивный алгоритм факторизации работает за время O(m^(1/2))
rho-алгоритм работает за время O(p^(1/2)), где p -- минимальный
простой делитель числа m => t <= O(m^(1/4)), поскольку p <= m^(1/2).
Пусть m -- 20-значное десятичное число. Тогда время работы <= 10^5.
=============================================
Колесо со спицами
Наивный алгоритм нахождения делителя числа m:
сначала делим m на 2
Если не делится, то делим на все пробные нечетные делители
3, 5, 7, 9, 11, 13, 15, 17, ...
до квадратного корня из m.
Если сначала поделим на 2 и на 3
То затем надо перебирать пробные делители, не кратные 2 и 3,
т.е. по модулю 6 они должны быть обратимы:
1, 5
Пробные делители
5,
7, 11,
13, 17,
19, 23,
25, 29,
. . .
Общая схема: берем набор небольших простых чисел, например,
smallPrimes = [2, 3, 5, 7, 11, 13]
их произведение равно
n = 30030
Сначала делим число m на 2, 3, 5, 7, 11, 13.
Если не делится ни на одно из них, то делим на следующие пробные
делители:
1) находим список чисел от 1 до n-1, взаимно простых c n:
inversesMod_n = [
d for d in range(1, n) if gcd(d, n) == 1
]
Длина этого списка равна 5760.
2) перебираем пробные делители d вида
d = k*n + s, (исключая d == 1),
где k = 0, 1, 2, 3, 4, ...
и для каждого фиксированного k числа s пробегают
все элементы списка inversesMod_n
Делители перебираются до квадратного корня
-------------------------------------------------------
p-1 алгоритм факторизации Полларда
Основан на Малой теореме Ферма.
p -- простое число, b != 0 (mod p) ==>
b^(p-1) == 1 (mod p)
Дано: число m
Надо: разложить m на множители.
Пусть p | m и число p-1 является гладким (smooth):
p - 1 = q1^e1 * q2^e2 * ... * q_k^e_k
s1 * s2 * ... * s_k
s_i = q_i^e_i
Пусть s_i <= N, где N -- параметр нашего алгоритма
например, N = 10,000,000
Возьмем случайное число b <- Z_m, проверим, что b взаимно просто с m
(если не так, то gcd(b, m) != 1 -- нетривиальный делитель m).
Тогда p | b^(N!), и поэтому p | gcd(b^(N!) - 1, m)
Доказательство:
p - 1 = s1*s2*...*s_k, s_i <= N, s_i -- степени разных простых чисел
Значит, s1*s2*...*s_k | N!
N! = (p-1)*t
По Малой теореме Ферма
b^(N!) = b^((p-1)*t) == (b^(p-1))^t == 1^t == 1 (mod p)
ПЕРЕРЫВ до 10:45
Замечание к p-1 алгоритму Полларда.
Мы возводили b в степень N!
Но можно на самом деле немного убыстрить алгоритм,
возводя не в степень N!, а в произведение всех чисел
s <= N,
которые являются максимальными степенями простых чисел
s = p^e <= N, s^(e+1) > N
Например, пусть N = 100
Тогда можно возводить в степень, равную произведению чисел
64, 81, 25, 79, 11, 13, 17, 19, ..., 97.
Также можно вычислять gcd не на каждом шаге, а, скажем, на каждом
10-м шаге (???).
===================================================
Серьезное усиление p-1 алгоритма Полларда было получено Ленстрами
с помощью использования эллиптических кривых над кольцом Z_m.
Эллиптическая кривая -- это кривая на плоскости R^2,
которая задается уравнением
y^2 = x^3 + a*x + b
Любая прямая пересекает элл. кривую в 3-х точках.
Обозначим их через p, q, r
Можно ввести операцию сложения точек элл. кривой:
p + q = r', где точка r' симметрична точке r относительно оси x:
r = (x, y)
r' = (x, -y)
Ноль -- это бесконечно удаленная точка кривой, в однородных координатах
(0, 1, 0) с x == 0
ВАЖНО!!!
Сумма точек на эллиптической кривой вычисляется с помощью четырех
арифметических операций +, -, *, /.
Относительно этой операции точки кривой превращается в абелеву группу.
Ноль -- это бесконечно удаленная точка кривой.
Малая теорема Ферма -- это теорема Лагранжа, примененная к группе обратимых
элементов кольца Z_p. Ее порядок равен p-1.
Порядок элемента делит порядок группы:
b^(p - 1) == 1 (mod p)
Давайте применим ту же теорему к группе точек эллиптической кривой
на полем Z_p. В случае элл. кривой мы используем аддитивную запись
В мультипликативной записи
b^k = b*b*b*...*b
k раз
В аддитивной записи
k*b = b+b+b+...+b
k раз
Операция умножения на k реализуется с помощью алгоритма,
аналогичного алгоритму быстрого возведения в степень.
Дано:
чило m, которое надо разложить на множители.
Выбираем случайную эллиптическую кривую и случайную начальную
точку b0 на ней.
И дальше вычисляем последовательность
b0,
b1 = 2*b0,
b2 = 3*b1,
b3 = 4*b2,
...
b_N = (N+1)*b_N-1
как и в p-1-алгоритме Полларда.
Пусть K = порядок группы точек эллиптической кривой над Z_p, где p | m.
Если K яляется N-гладким числом, то b_N равняется нулевой точке элл. кривой
=> координата x этой точки равну нуля (mod p) => координата делится
на p и dgc(x, m) делится на p.
Зачем всё это нужно?
Представим, что p-1 не является гладким числом => m нельзя разложить
с помощью p-1-алгоритма (нам не повезло).
Но с эллиптическими кривыми число точек элл. кривой зависит от чисел a, b,
порядок группы точек не постоянен. Если нам не повезло с какой-то кривой,
то выбираем другую (случайную) кривую.
Реализация алгоритма:
при вычислении суммы точек мы используем деление в Z_m. Но если элемент
x необратим в Z_m, то x не взаимно прост с m,
gcd(x, m) != 1
и мы разложили число m на множители.
Для нулевой точки p = (x, y, z) x == 0 (mod p)
и p | gcd(x, m) => gcd(x, m) -- нетривиальный делитель m.
=========================================================
Алгоритм вычисления квадратного корня в поле вычетов по модулю p
(p -- простое число).
Дано: p -- простое число
a -- элемент поля Z_p, являющийся квадратом (что легко проверить)
Надо: вычислить квадратный корень x из a:
x^2 == a (mod p)
Предложение:
элемент a != 0 является квадратом (mod p) <==>
a^((p-1)/2) == 1 (mod p)
Доказательство вытекает из того, что группа обратимых элементов любого поля
является циклической.
Пусть x -- притивный корень порядка p-1 (порождающий элемент
группы ненулевых элементов поля Z_p).
Для любого a != 0 (mod p)
a = x^e (mod p)
a есть квадрат <=> e четное число
Пусть a -- квадрат, т.е. a = x^e, e -- четное.
a^((p-1)/2) = (x^e)^((p-1)/2) = (x^(e/2))^(p-1) == 1 (mod p)
Если a -- не квадрат, то e нечетное и
в циклической группе элементы нечетной степени отличны от 1.
a^((p-1)/2) = (x^(2s + 1))^((p-1)/2) =
x^((2*s + 1)*((p-1)/2)) == x^((p-1)/2) != 1
поскольку x порождает группу ненулевых элементов Z_p
и его порядок в точности равен p-1.
(2*s + 1)*((p-1)/2) = 2s*(p-1)/2 + (p-1)/2
-------------------------------------------
Считаем, что p != 2 ==> p нечетное число.
Тогда p == 1 (mod 4) или p == -1 (mod 4) (или, что то же самое, p == 3 (mod 4))
Случай p == -1 (mod 4) совсем простой. Можно сразу
предъявить решение:
x == a^((p+1)/4) (mod p)
Проверим:
x^2 = a^((p+1)/2) = a^((p-1)/2)*a
поскольку (p+1)/2 = ((p-1) + 2)/2 = (p-1)/2 + 1
Так как a -- квадрат, то a^((p-1)/2) == 1 (mod p) и
x^2 = a^((p+1)/2) = a^((p-1)/2)*a == a (mod p)
Сложный случай: p == 1 (mod 4)
Число (p+1)/4 уже не является целым, и мы не может возвести a
в эту степень.
(p-1)/2 -- четное число.
Представим его в виде
(p-1)/2 = q * 2^s, где q -- нечетное число.
Затем вычисляем последовательность (mod p):
x0 = a^q,
x1 = x0^2,
x2 = x1^2,
...
x_s = x_s-1^2 = a^((p-1)/2) = 1 (mod p)
Последовательность имеет один из двух видов:
1, 1, 1, ..., 1
*, *, ..., *, -1, 1, ..., 1
Первый случай опять совсем простой:
a^q == 1, где q -- нечетное число
Тогда
r = a^((q + 1)/2)
является искомым квадратным корнем:
r^2 = a^(q + 1) = a^q * a = a (mod p),
задача решена.
Осталось рассмотреть последний случай, он самый сложный:
a^q != 1, последоватльеность имеет вид
*, *, ..., *, -1, 1, ..., 1
Обозначимs
t = a^q
r = a^((q + 1)/2)
Получаем
r^2 = a^(q + 1) = a^q * a = t*a
Инвариант цикла:
t -- корень степени 2^k из 1
Последовательность t, t^2, t^4, t^8, ..., 1
заканчивается единицей
Длина этой последовательности <= s
r^2 == t*a (mod p)
Цель: подправить t и r так, чтобы длина последовательности
t, t^2, t^4, t^8, ..., 1 (до единицы)
уменьшилась.
Окончание алгоритма будет в следующий раз.
---------------------------------------
18 Mar 2022
В прошлый раз мы рассматривали алгоритм извлечения квадратного
корня в поле вычетов по модулю p (p -- большое простое число).
Предложение. Пусть p != 2. Элемент a <- Z_p является квадратом <==>
a^((p-1)/2) == 1 (mod p)
Обычно такие числа называют "квадратичными вычетами по модулю p",
а элементы, которые не являются квадратами
z: z^((p-1)/2) == -1 (mod p)
назыают "квадратичными невычетами".
Утверждение:
группа ненулевых элементов конечного поля по умножению является циклической.
Пусть a -- ее порождающий элемент. Тогда любой x <- Z_p, x != 0
является некоторой степенью элемента a:
x = a^k
Элемент x является квадратом <==> k четное.
Число квадратичных вычетов, отличных от 0, равно числу квадратичных невычетов.
=> найти какой-нибудь квадратичный вычет или невычет не представляет
трудности (при поиске невычета обычно перебираем все элементы,
начиная с 2).
Алгоритм нахождения квадратного корня из a по модулю p.
Дано: p -- простое число, p != 2
a -- квадратичный вычет, a^((p-1)/2) == 1 (mod p)
Надо: найти r: r^2 == a (mod p)
Рассматриваем 2 случая
1. p == 3 (mod 4) (или, что то же самое, p == -1 (mod 4)).
Тогда мы сразу можем предъявить решение:
r == a^((p + 1)/4) (mod p)
Действительно,
r^2 == a^((p+1)/2) = a^((p-1)/2 + 1) = a^((p-1)/2) * a == a (mod p)
(p+1)/2 = ((p-1) + 2)/2 = (p-1)/2 + 1
Пример: p = 945152109870139
p == 3 (mod 4)
a = 234945000550770 -- это квадратичный вычет,
a = 1234567891011^2 (mod p)
Вычисляем r:
r == a^((p+1)/4) = 943917541979128
-r == p - 234945000550770 == 1234567891011 (mod p)
Остается рассмотреть сложный случай: p == 1 (mod 4)
В этом случае p-1 делится на 4, значит, (p-1)/2 -- четное число.
Выделим в нем максимальную степень двойки:
(p - 1)/2 = q*2^s, где q -- нечетное число
Рассмотрим последовательность длины s+1:
t == a^q (mod p)
t^2, (mod p)
t^4, (mod p)
t^8,
...
t^(2^s) = 1 (mod p)
В этой последовательности каждый следующий элемент есть квадрат предыдущего.
Она имеет один из двух видов: если t == a^q == 1 (mod p), то
1, 1, 1, ..., 1 (длина s+1)
или, когда t == a^q != 1 (mod p),
*, *, ..., *, -1, 1, ..., 1
В первом случае имеем
t == a^q == 1 (mod p), q -- нечетное.
Положим
r == a^((q+1)/2)
Это и будет искомый корень!
r^2 == a^(q+1) == a^q * a == 1 * a == a (mod p)
И теперь самый сложный случай:
a^q != 1 (mod p)
Мы обозначили
t == a^q (mod p)
Точно так же вычислим r:
r == a^((q+1)/2) (mod p)
Тогда r^2 равно
r^2 == a^(q+1) == a^q * a == t*a (mod p)
Итак, сформулируем инвариант цикла для пары (t, r):
t -- это корень из 1 степени 2^s,
последовательность квадратов, которая начинается с t
и имеет длину s+1, заканчивается единицей
t, *, ..., *, -1, 1, ..., 1
r^2 == t*a (mod p)
Идея: подправим эти два числа t, r так, чтобы инвариант сохранился,
но число единиц в конце последовательности увеличилось. В конце-концов
на некотором шаге t станет равным единице.
Как подправить эти два элемента?
Здесь мы используем любой квадратичный невычет z:
z не является квадратом (mod p)
z^q, z^(2*q), z^(4*q), ..., z^((p-1)/2) == -1 (mod p)
Последовательность квадратов для z^q имеет вид
*, *, ..., *, -1 (длина последовательности s+1)
Последовательность квадратов для t:
*, *, ..., *, -1, 1, ..., 1, 1
Для z^q
*, *, ..., *, *, *, ..., *, -1
Имеем
t
*, *, ..., *, -1, 1, ..., 1, 1
*, *, ..., *, *, *, ..., *, -1
b, b^2
Модифицируем t и r следующим образом:
t' == t*b^2 (mod p)
r' == r*b
Тогда инвариант r^2 == t*a (mod p) сохранится:
r'^2 = (r*b)^2 == r^2 * b^2 == t*a * b^2 == t'*a
Но последовательность квадратов для t' до -1 имеет меньшую длину,
поскольку она равна произведению последовательность квадратов для t
на последовательность квадратов для b^2, где два элемента -1
перемножаются и дают 1. Значит, не более чем за s шагов t станет равным
единице и r будет искомым квадратным корнем
(т.к. r^2 == t*a (mod p), t == 1 ==> r^2 == a).
Рассмотрим на примере:
p = 7463806529
(p - 1)/2 = q * 2^s, где
q = 116621977, s = 5
Пусть a = 670548136 (это квадратичный вычет)
t == a^q == 1450582071
r == a^((q+1)/2) == 711481791
Проверим инвариант r^2 == t*a:
r^2 == 868383522
t*a == 868383522 всё верно.
Находим квадратичный невычет z = 3
Последовательность квадратов для x == z^q (mod p):
7352389316, 7261267627, 627972008, 5676044929, 2110798942, -1, 1
Последовательность квадратов для t (mod p):
1450582071, 1787761600, 2110798942, -1, 1
Корректирующий множитель b = 7261267627, b^2 = 627972008
Корректируем пару (t, r):
t' == t*b^2 (mod p) == 3183782943
r' == r*b (mod p) == 1618420748
Выписываем последовательность квадратов для t':
3183782943, 5353007587, -1, 1
Находим корректирующий множитель:
7352389316, 7261267627, 627972008, 5676044929, 2110798942, -1, 1
=========
b b^2
b = 627972008
t'' = t'*b^2 (mod p) == 1
r'' = r*b (mod p) == 7340349740
Получили ответ: r'' == 7340349740 (mod p)
Проверка:
(r'')^2 == 670548136 (mod p)
a == 670548136 (mod p)
Задача решена правильно.
ПЕРЕРЫВ до 10:50
Классы в языке Python
Объектно-ориентированные языки
SmallTalk (фирма Xerox)
C++ (не является объектно-ориентированным языком в строгом смысле,
занимает промежуточное положение между традиционными языками,
такими как Fortran, C, Pascal, и настоящими
объектно-ориентированными языками, такими как SmallTalk, Forte,
Java, C# (и другие языки платформы .Net, Visual Basic, CLR C++),
Python, ...)
Java (начало 90-х годов, первоначальное название Oak -- язык задумывался как
объектно-ориентированный вариант C++; в Java sin -- это статический
метод класса Math, т.е. надо писать Math.sin(x)).
C# (.Net)
. . .
Python
Черта, специфическая для объектно-ориентированных языков -- это наличие
объектов. Все переменные устроены одинокаво (пусть даже у них может быть
тип, как в Java/С#), они не хранят объекты, как в С++, а хранят
объектные ссылки. Все объекты живут в контролируемой динамической памяти.
Объекты нельзя удалять самому, этим занимается специальный процесс
gc -- garbage collector. Он МОЖЕТ удалить только беспризорные объекты,
на которые уже нет ссылок и которые никому не нужны.
Например, вот такой код на C++ -- это ГРУБАЯ ОШИБКА:
A& f() {
A res;
. . .
return res;
}
Возвращается ссылка на локальный объект, который удаляется
при выходе из функции. ТАК ПИСАТЬ НЕЛЬЗЯ!!!
А в языке Java так писать можно:
A f() {
A res = new A();
. . .
return res;
}
x = f();
Переменные в объектно-ориентированных языках содержат
объектные ссылки, а не сами объекты, как в С++.
Особенность Питона по сравнению с большинством других
ООЯ: Питон -- это язык с динамической типизацией.
Переменные не имеют типов, вообще в Питоне нет описаний
типов! (Очень непривычно, потому что всегда раньше говорилось,
что уровень языка тем выше, чем больший процент в нем занимают
описания.)
В С++ описание класса содержит как бы чертеж объекта.
Все объекты в C++ имеют одинаковое устройство.
В Питоне это СОВЕРШЕННО НЕ ТАК!
Устройство двух объектов одного и того же класса может быть
разным и зависеть от истории их создания и модификации
(конечно, программисты стараются так не делать, но формально
ограничений нет).
Переменные-члены класса в Питоне создаются при вызове методов класса,
а не описываются в h-файле, как в C++.
Рассмотрим устройство класса на примере класса "Вектор на плоскости"
R2Vector.
=============================================================
25 Mar 2022
Классы в Питоне
Классы в Питоне не похожи на классы в C++.
Главное отличие проистекает из того, что Python -- язык с динамической
типизацией. Если в C++ мы определяем как бы "чертеж" класса (строение
его объектов)
class R2Vector {
public:
double x;
double y;
public:
R2Vector(double xx = 0., double yy = 0.):
x(xx),
y(yy)
{}
. . .
double norm() const {
return sqrt(x*x + y*y);
}
/* Так писать не надо!
R2Vector operator+(const R2Vector& v) const {
return R2Vector(x + v.x, y + v.y);
}
*/
R2Vector& operator+=(const R2Vector& v) {
x += v.x;
y += v.y;
return *this;
}
// Умножение: 2 варианта
// Это либо умножение вектора на число,
// либо скалярное произведение двух векторов.
/* Вместо этого метода используем глобальную функцию
заданную вне класса
R2Vector operator*(double c) const {
return R2Vector(x*c, y*c);
}
*/
double operator*(const R2Vector& v) const {
return x*v.x + y*v.y;
}
R2Vector& operator*=(double c) {
x *= c;
y *= c;
return *this;
}
. . .
};
inline R2Vector operator+(const R2Vector& u, const R2Vector& v) {
R2Vector res(u);
res += v;
return res;
}
Плохой дизайн языка С++: все время приходится реализовывать смесь
классов и глобальных функций. Невозможно обойтись только классами.
Например, в языке Java вообще нет функций!
y = sin(x); // Так написать в Java нельзя!
y = Math.sin(x); // Надо писать так.
// Здесь Math.sin(x) -- это статический метод
// класса Math
Мы написали з метода умножения вектора как методы класса.
Ими можно пользоваться так:
R2Vector u, v;
double c;
double s = u*v;
u *= c;
R2Vector w = 2.5*u; // Ошибка! Первый аргумент вообще не является классом
// Надо писать так:
R2Vector w = u*2.5; // Коэффициент обязательно должен быть
// правым аргументом операции умножения *
y = f(x) левая запись отображения
y = (x)f правая запись отображения -- соответствует
объектно-ориентируемому стилю.
В С++ приходится вместо метода класса "operator*"
использовать две глобальных функции
R2Vector operator*(double c, const R2Vector& v) {
return R2Vector(c*v.x, c*v.y);
}
R2Vector operator*(const R2Vector& v, double c) {
return R2Vector(v.x*с, v.y*с);
}
В Питоне можно внутри класса реализовать как правое умножение
вектора на число, так и левое умножение:
class R2Vector:
def __init__(self, x = 0., y = 0.):
if type(x) == list or type(x) == tuple:
self.x = float(x[0]);
self.y = float(x[1]);
else:
self.x = float(x);
self.y = float(y);
def __mul__(self, c): # Умножение вектора self СПРАВА
# на вектор c или число c
# w = v*c
if type(c) == R2Vector or type(c) == list or type(c) == tuple:
return self.x*float(c[0]) + self.y*float(c[1])
else:
return R2Vector(self.x*float(c), self.y*float(c))
def __rmul__(self, c): # Умножение вектора self СЛЕВА
# на вектор c или число c
# w = c*v
if type(c) == R2Vector or type(c) == list or type(c) == tuple:
return float(c[0])*self.x + float(c[1])*self.y
else:
return R2Vector(self.x*float(c), self.y*float(c))
==============================================================
Переменные класса и статические методы класса
Пример: мы хотим реализовать кольцо вычетов по модулю m
class Zm:
m = 7
@staticmethod
def setmod(mod):
Zm.m = mod
@staticmethod
def getmod(mod):
return Zm.m
__init__(self, x = 0):
self.n = x%Zm.m
======================================================
Перерыв 10 минут до 10:35
Работа с матрицами в Питоне. Модуль numpy, работа с матрицами и
многомерными массивами в numpy
Питон -- это скриптовый язык: он очень медленный и серьезные
вычисления в нем не делают; вместо этого используются модули Питона,
которые на самом деле написаны на традиционных языках программирования:
Fortran, C, C++. В результате на Питоне можно писать настоящие программы,
а не только учебные или игрушечные. Например, для умножения матриц A и B
в самом
Питоне мы исполняем только одну команду
C = numpy.dot(A, B)
а само умножение уже выполняется быстро внутри модуля numpy,
который написан уже не на Питоне.
Кстати, умножение матриц (или скалярное произведение векторов)
в Питоне можно записывать с помощью операции @
С = A@B
План:
реализуем алгоритм Гаусса приведения матрицы к ступенчатому виду
(row echelon form, ref) для матриц, определенных в стиле Питона,
с вещественными элементами.
Алгоритм гаусса приведения матрицы к ступенчатому виду
Используются элементарные преобразования строк матрицы 2-х типов:
1) меняем местами две строки матрицы с индексамим i, k;
это преобразование меняет знак определителя. Если мы хотим сохранить
определитель, то можно при этом менять знак одной из строк;
2) к i-й строке прибавляем k-ю строку, умноженную на коэффициент c;
3) умножить i-ю строку матрицы на коэффициент c.
Инвариант цикла в алгоритме Гаусса:
задается двумя индексами i, j
Минор матрица в первых j столбцах уже приведен к ступенчатому
виду, его ранг (число ненулевых строк) равен i.
Необработанная часть матрицы представляет собой минор,
начинающийся со строки i и столбца j
1. Ищем максимальный по модулю элемент в j-м столбце, начиная со строки i
. . .
С помощью метода Гаусса можно
1) решить систему линейных уравнений;
2) вычислить обратную матрицу.
1. Решение системы с невырожденной матрицей.
AX = B
В -- вектор-столбец свободных членов.
Мы выписываем расширенную матрицу системы
[1 2 3] [2]
A = [4 5 6] B = [3]
[7 8 0] [5]
Формируем расширенную матрицу:
[1 2 3 | 2]
[4 5 6 | 3]
[7 8 0 | 5]
Приводим расш. матрицу к ступенчатому виду:
[* * * | *]
[0 * * | *]
[0 0 * | *]
Выполняем обратный проход:
1) идем от последней к первой строке
2) очередную строку домножаем на элемент, обратный к диагональному;
[* * * | *]
[0 * * | *]
[0 0 1 | *]
3) из вышестоящих строк с индексом k вычитаем текущую строку i, умноженную на
коэффициент a_ki
[* * 0 | *]
[0 * 0 | *]
[0 0 1 | *]
В результате зануляется столбец i выше очередного диагонального элемента.
В результате обратного прохода получаем матрицу
[1 0 0 | *]
[0 1 0 | *]
[0 0 1 | *]
Столбец свободных членов и будет решением системы.
Аналогично вычисляется обратная матрица
Формируем расширенную матрицу:
[1 2 3 | 1 0 0]
[4 5 6 | 0 1 0]
[7 8 0 | 0 0 1]
Приводим к ступенчатому виду:
[* * * | * * *]
[0 * * | * * *]
[0 0 * | * * *]
Выполняем обратный проход:
[1 0 0 | * * *]
[0 1 0 | * * *]
[0 0 1 | * * *]
Справа будет обратная матрица.
Замечание:
нельзя на практике вычислять обратную с помощью алгебраических дополнений!!!
A_ij (алгебраическое дополнение) -- определитель минора, который
получается вычеркиванием из матрицы a i-й строки и
j-го столбца,
Элементы обратной матрицы равны (-1)^(i+j) A_ji/det(a)
По этой формуле в принципе можно вычислить обратную матрицу. Но
1) сколько операций потребуется??? O(n^5)
A_ij вычисляем за O(n^3)
всего n^2 элементов, получается n^2 * O(n^3) = O(n^5)
Для матрицы порядка 100 алгоритм будет работать в 10000 раз дольше,
чем рассмотренный выше алгоритм, основанный на методе Гаусса и
обратном проходе (заботает за O(n^3));
2) ничего не понятно с устойчивостью этого метода.
Домашнее задание: выполнить подобные задачи для
матриц, представленных как numpy-array.
====================================
1 Apr 2022
Работа с матрицами в Питоне
Прошлый раз мы представляли матрицу, используя базовый тип
"список" (list) Питона. Матрица представлялась списком строк,
каждая строка -- список чисел
>>> a = [[1., 2., 3.], [4., 5., 6.], [7., 8., 0.]]
>>> x = a[1][1]
>>> x == 5.
True
В реальных задача никто так матрицы не представляет!
Почему? Программа будет работать очень медленно, в тысячи или миллионы
раз медленнее, чем такая же программа на C++.
Надо понимать, что Python -- это "скриптовый" язык, т.е.
мы на Питоне пишем только программы, состоящие из команд самого
верхнего уровня. При этом мы всегда используем модули ("батарейки"),
которые работают с разнообразными объектами, написанные на "настоящих"
яызках программирования: Fortran, С, С++. Напимер, матрица представляется
как элемент класса numpy.ndarray модуля numpy, который, конечно же,
написан не на Питоне и в котором все действия с матрицами реализованы
очень эффективно.
import numpy as np
a = np.array(...)
b = np.array(...)
c = np.dot(a, b)
c = a @ b # def __matmul__(self, b):
Для массивор из модуля numpy арифметические операции определны как
векторные, это сильно отличается от операций со списками Питона.
Операции, как правило, выполняются покомпонентно.
Например, умножение двух массивов одинаковой длины выполняется покомпонентно,
результатом является массив той же длины. Для скалярного произведения
векторов, произведения матриц и т.п. имеется функция
c = np.dot(a, b)
(скалярное произведение == dot product).
Ее в современном Питоне можно обозначать операцией @
c = a @ b
Работа с матрицами
Матрица представляется как объект типа numpy.ndarray,
в котором хранится матрица в стиле языка C:
1) все элементы матрицы имеют одинаковый тип;
2) они хранятся в непрерывном сегменте памяти.
Достаточно необычно реализован доступ к элементам матрицы.
Он позволяет выполнять векторные операции со строками, столбцами
матрицы, используя одну команду Питона.
>>> import numpy as np
>>> a = np.array([[1., 2., 3.], [4., 5., 6.], [7., 8., 0.]])
>>> a
array([[1., 2., 3.],
[4., 5., 6.],
[7., 8., 0.]])
>>> a.shape
(3, 3)
Строка матрицы с индексом 1:
>>> a[1]
array([4., 5., 6.])
Диагональный элемент матрицы с индексами 1, 1
>>> a[1, 1]
5.0
Минор матрицы, состоящий из строк с индексами 0, 1
>>> a[[0, 1]]
array([[1., 2., 3.],
[4., 5., 6.]])
Обмен двух строк матрицы с индексами 0, 1
>>> a[[0, 1]] = a[[1, 0]]
>>> a
array([[4., 5., 6.],
[1., 2., 3.],
[7., 8., 0.]])
К строке матрицы с индексом 0 прибавить строку с индексом 1,
умноженную на число 10.:
>>> a
array([[4., 5., 6.],
[1., 2., 3.],
[7., 8., 0.]])
>>> a[0] += a[1]*10.
>>> a
array([[14., 25., 36.],
[ 1., 2., 3.],
[ 7., 8., 0.]])
Умножить на 10. строку матрицы с индексом 2:
>>> a
array([[14., 25., 36.],
[ 1., 2., 3.],
[ 7., 8., 0.]])
>>> a[2] *= 10.
>>> a
array([[14., 25., 36.],
[ 1., 2., 3.],
[70., 80., 0.]])
ПЕРЕРЫВ 10 минут до 10:35
Файл npmatrix.py содержит реализацию метода Гаусса
для матриц, представленных как объекты типа numpy.ndarray
(так принято представлять матрицы во всех реальных программах
на Питоне).
==========================================
Оконная графика на Python: модуль tkinter
tkinter -- простейший модуль для оконного программирования.
Его можно быстро выучить и писать несложные оконные программы.
==============================================
8 Apr 2022
Оконное и графическое программирование на Python:
модуль tkinter
Название произошло от пакета TCL/TK, использовался в Unix-системах.
Я в начале использовал модуль kde (Кommon Desktop Environment)
=
X-Window -- распределенная оконная система, где на X-терминалах
работает процесс XServer, а на главном компьютере запускаются
оконные приложения, а также программа, которая называется
Window Manager. Она отвечает за расположение окон на экране X-терминала,
иконки на поверхности рабочего стола, toolbar'ы, taskbar и т.п.
Несколько фирм создали общий оконный интерфейс (i.e. Window Manager)
под названием CDE (Common Desktop Environment) -- но это не был открытый
продукт, а идеология GNU (Gnu Is Not Unix).
Все коммерческие утилиты Unix были переписаны как свободные варианты,
которые, как правило, намного лучше коммерческих. Названия ироничные:
more (просмотр файлов) -- less
yacc (Yet Another Compiler Compiler инструмент для создания парсеров) --
bison
cde -- kde
OK -- All Correct
Сначала я не знал, что в базовом дистрибутиве MS Windows есть
модуль для разработки оконных приложений (я пользовался модулем kde,
которого нет в базовой установке Python под MS Windows).
Но оказалось, что там всё-таки есть модуль tkinter, который позволяет
быстро и легко создавать (не самые сложные) оконные приложения.
Плюсы tkinter:
1) очень простой оконный интерфейс, буквально сделан для детей;
2) совершенно независим от операционной системы, компьютера и т.д.
Недостатки:
1) под X-Window не поддерживает сглаживание линий
(antialiasing), true-type fonts (???) и т.п., только базовые
возможности X-Windows. Под Mac OS/X таких проблем нет, поскольку
там используется Quartz, в котором всё это реализовано уже на уровне
операционной системы;
2) возможно, в tkinter'е нет каких-нибудь продвинутых оконных
интерфейсов.
Если нужно что-то быстро написать (оконный вариант программы),
то tkinter -- хороший выбор.
ПЕРЕРЫВ 12 минут до 10:45
. . .
==========================================
15 Apr 2022
Метод наименьших квадратов
Задача регрессии
Пусть надо построить зависимость зависимой величины y <- R
от независимых величин x1, x2, ..., x_n <- R. Числа x_i -- это
признаки некорого объекта (features), y = y(x1, x2, ..., x_n) --
некоторая характеристика объекта. Мы хотим построить
линейную модель, которая вычисляет зависимую величину y
по независимым величинам x1, ..., x_n. Рассмотрим этот набор
как вектор из R^n:
x = (x1, x2, ..., x_n) <- R^n
Мы ищем линейную зависимость
y = f(x) = <w, x> + b, w = (w1, w2, ..., w_n) <- R^n, b <- R.
Угловыми скобками обозначено скалярное произведение векторов.
Пусть мы имеем тренировочную выборку из m объектов
x1 = (x_11, x_12, ..., x_1n), y1 = y(x1),
x2 = (x_21, x_22, ..., x_2n), y2 = y(x2),
. . .
x_m = (x_m1, x_m2, ..., x_mn), ym = y(xm).
Задача состоит в том, чтобы построить линейную зависимость по
этой выборке, т.е. найти параметры w, b, которые обеспечивают
зависимость, максимально близкую к этой выборке.
f(x) = f_w,b(x) = <w, x> + b
Заметим, что можно добавить еще один признак 1 ко всем объектам
выборки
x_i' = (x_i1, x_i2, ..., x_in, 1) <- R^(n+1)
Обозначим w' = (w, w0), где w0 = b
Тогда f(x) выражается в виде скалярного произведения:
f(x') = <w', x'>
Задача упрощается.
Давайте опустим штрих,
Мы строим модель, которая минимизирует среднюю ошибку.
В методе наименьших мы минимизируем среднюю квадратичную ошибку
MSE -- Mean Squared Error:
MSE(w) = 1/m Sum_i (f(x_i) - y_i)^2 --> min_w
MSE(w) = 1/m Sum_i (<w, x_i> - y_i)^2 --> min_w
Запишем эту задачу в матричном виде:
[x_11, x_12, ..., x_1n] [ y1 ]
A = [x_21, x_22, ..., x_2n] Y = [ y2 ] m >= n
[. . . ] [... ]
[x_m1, x_m2, ..., x_mn] [ y_m]
[ w1 ]
W = [ w2 ]
[ ... ]
[ w_n ]
В матричном виде мы как бы решаем уравнение
AW = Y
где A -- это прямоугольная матрица размера mxn,
W -- вектор-столбец размера n,
Y -- вектор-столбец размера m.
Мы получаем систему, в которой число уравнений m больше, чем
число неизвестных n. Она в общем случае не имеет решения,
но мы как "заменитель" решения ищем вектор V такой, что
квадрат евклидова расстояния между векторами AV и Y
был бы минимальным.
A -- это линейный оператор
A: R^n --> R^m, m >= n
Образ этого линейного оператора -- это подпространство Im(A) < R^m
Надо найти вектор V <- R^n такой, что Y* = AV лежит максимально
близко к вектору Y
т.е. надо минимизировать расстояние между векторами Y* = AV и Y.
Понятно, что Y* равняется проекции вектора Y на подпространство
Im(A), т.е. разность векторов Y - Y* перпендикулярна
подпространству Im(A), т.е. перпендикулярна любому вектору из Im(A)
<Y - AV, Aw> = 0 для любого вектора w <- R^n
<Y, Aw> - <AV, Aw> = 0
Y^t A w = (AV)^t A w
Поскольку это матричное равенство выполняется для любого вектора w,
то матрицы, которые умножаются на w, равны:
Y^t A = (AV)^t A
Y^t A = V^t A^t A
Транспонируем это равенство:
A^t Y = A^t A V
В линейной алгебре хорошо известно утверждение:
если столбцы прямоугольной матрицы A линейно независимы,
то матрица A^t A обратима (это квадратная матрица размера nxn).
Домножим равенство на обратную к ней матрицу (A^t A)^(-1)
(A^t A)^(-1) A^t Y = V
Мы нашли вектор V.
Эта матрица
A+ = pinv(A) = inv(A*A) A* (звездочкой здесь
мы обозначили транспонирование)
называется псевдообратной матрицей Мура--Пенроуза,
она обощает понятие обратной матрицы для прямоугольной матрицы A.
В общем случае псевдообратная матрица вычисляется с помощью
сингулярного разложения матрицы:
A -- любая прямоугольная матрица. Она представляется
в виде
A = U D V,
где U, V -- ортогональные квадратные матрицы UU* = 1, VV* = 1
D -- диагональная прямоугольная матрица.
(Это называется SVD-разложение, Singular Value Decomposition).
Тогда в общем случае псевдообратная матрица A+ выражается как
A+ = V* inv(D) U*
В Питоне все эти алгоритмы реализованы в модуле numpy.linalg
(подмодуль модуля numpy, название означает "линейная алгебра").
ПЕРЕРЫВ 10 минут до 10:40
Давайте реализуем оконную программу на Питоне:
мы кликами мыши отмечаем узлы интерполяции, а также задаем степень
полинома с помощью слайдера. По нажатию на кнопку Draw
программа вычисляет коэффициенты полинома, используя метод наименьших
квадратов, и рисует график полинома, минимально уклоняющегося
от заданных узлов.
Используем в качестве образца программу plot_func.py,
которую мы написали на прошлом занятии.
Программа leastSquares.py
. . .
Вопрос: почему задача нахождения зависимости между независимой и
зависимой величинами называется задачей регрессии?
Был ученый, который попытался определить зависимость
между ростом отца и ростом его (взрослого) сына.
Получаем набор точек p, где p.x -- рост отца, p.y -- рост сына.
После построения прямой по набору точек получилось, что угол наклона
прямой меньше 45 градусов. Это означает, что, если отец высокий,
то сын в среднем ниже его и, обратно, если отец низкий, ещ сын
в среднем его выше. То есть человеческое общество в среднем
эволюционирует в направлении усреднения роста (деградация,
регрессия).
Еще: при машинном обучении мы всегда строим модель, минимизируя
функцию ошибки (Error Function, Loss Function) на некоторой
тренировочной выборке. В методе наименьших квадратов
мы минимизируем сумму квадратов отклонений реальных измерений
от значений, которые выдает наша модель:
MSE(w) = 1/m Sum_i (f(x_i) - y_i)^2 --> min_w
Здесь f(x) -- это значение, которая наша модель выдает на
объекте x, а y_i -- реальное значение, измеренное на
конкретном объекте x_i
В качестве меры ошибки (Loss Function)на одном объекте
используется квадрат разности:
L(y', yy'') = |y' - y''|^2
Эта функция плохая, потому что она очень чувствительна
к выбросам (outliers).
Если вместо квадрата разности испольпользовать пролсто абсолютную
величину разности
E(y', yy'') = |y' - y''|
то выбросы будут влиять на ответ гораздо меньше (все точки имеют
одинаковое на величину ошибки).
Такая функция ошибки называется MAE (Mean Average Error):
MAE(w) = 1/m Sum_i |f(x_i) - y_i| --> min_w
В чем проблема? Функция y = |x| не дифференцируема в нуле,
а точного решения, такого как в методе наименьших квадратов,
при минимизации функции MAE не существует. Поэтому надо искать минимум
функции ошибки MAE, используя какие-нибудь итеративные методы
минимизации: градиентный спуск, метод сопряженных градиентов
(вернее, итеративный метод на его основе), а все эти методу
хорошо работают только для дифференцируемых функций.
Выход простой: используем функцию Хубера (Huber's function)
{ 1/2 * x^2 при |x| <= delta
H(x) = {
{ delta*(|x| - delta/2) при |x| > delta
Модель строится путем минимизации функции ошибки, использующей
функцию Хубера:
HuberLoss(w) = 1/m Sum_i H(f(x_i) - y_i) --> min_w
В Питоне минимизация функции ошибки, использующей функцию Хубера,
реализована в пакете sklearn:
import numpy as np
from sklearn.linear_model import HuberRegressor
def onDraw():
. . .
a = np.vander([points[i].x for i in range(m)], degree + 1)
y = np.array([points[i].y for i in range(m)])
huber_reg = HuberRegressor(fit_intercept=False)
huber_reg.fit(a, y)
w = [huber_reg.coef_[i] for i in range(degree + 1)]
Домашнее задание
1. Реализовать оконную программу, по данному набору узлов
и заданной степени многочлена вычислить многочлен методом
наименьших квадратов (минимизирующий функцию ошибки MSE)
и нарисовать его график.
2.Реализовать оконную программу, по данному набору узлов
и заданной степени многочлена вычислить многочлен путем
минимизации функции потерь Хубера (приблизительно
минимизируя функцию ошибки MAE) и нарисовать его график.
3. Реализовать оконную программу, по данному набору узлов
и заданной степени многочлена вычислить многочлен методом
линейного программирования, минимизирующий функцию ошибки MAE,
и нарисовать его график.
=======================================
22 Apr 2022
Была рассмотрена задача линейной регрессии
У нас есть независимая величина x <- R^n и зависящая от нее
величина y. У нас есть набор измерений величины y для различных
значений x:
x1, x2, x3, ..., x_m
y1, y2, y3, ..., y_m
Откуда взялось слово "регрессия"?
Ученый изучал зависимость роста (взрослого) сына от роста отца.
Для конкретной выборки от изобразил пару (рост отца, рост сына)
точкой на плоскости: по оси X рост отца, по Y -- сына.
Дальши он провел прямую, наиболее близкую к этим точкам,
и получилось:
1) угол наклона прямой меньше 45 градусов;
2) прямая проходит не через начало координат, а немного выше,
через точку (0, y0), где y0 > 0.
Это означает, что у высокого отца сын чаще всего ниже его,
а у низкого отца -- наоборот, выше.
Вывод: со сменой поколений рост усредняется, т.е. происходит
"регрессия".
Линейная модель:
f(x) = <w, x> + b
Например, как строится полином степени d в методе наименьших кадратов?
Есть точки z1, z2, ..., z_m <- R
По ним строим наборы степеней
x1 = (z1^d, z1^d-1, ..., z1, 1)
x2 = (z2^d, z2^d-1, ..., z2, 1)
x3 = (z3^d, z3^d-1, ..., z3, 1)
. . .
x_m = (z_m^d, z_m^d-1, ..., z_m, 1)
Пусть многочлен имеет коэффициенты w = (w0, w1, ..., w_d)
p(z) = w0*z^d + w1*z^d-1 + ... + w_d-1*z + w_d =
= <w, x>, где x = (z^d, z^d-1, ..., z, 1)
В данном случае свободный член b = 0.
От свободного члена всегда можно избавится, добавив к описанию объекта
дополнительный произнак, тождественно равный 1.
x1, x2, ..., x_m <- R^n
(x1,1), (x2,1), ..., (x_m,1) <- R^n+1
Для исходных объектов мы имели зависимость со свободным членом
(intercept) b:
f(x) = <w, x> + b
Для "расширенных" объектов тогда будет зависимость без свободного
члена:
g(x') = <w', x'>,
где x' = (x, 1), w' = (w, b)
В прошлый раз было рассказано:
1) метод наименьших квадратов имеет точное решение:
ищем зависимость
f(x) = <w, x>
которая минимизирует средне-квадратичную ошибку
MSE = 1/m Sum_i (f(x_i) - y_i)^2
В векторном виде мы как бы ищем решение системы линейных
уравнений
A W = Y
где A -- это матрица, построенная по векторам x_i, i = 1,...,m
x_11, x_12, ..., x_1n
x_21, x_22, ..., x_2n
. . .
x_m1, x_m2, ..., x_m,n
причем m >= n
Y -- столбец
y1
y2
...
ym
Число уравнений m больше (или равно), чем число неизвестных =>
точного решения системы нет.
A W = Y
Но мы ищем такое решение V, что разность между левой и правой
частями уравнения была бы минимальной:
A V = Y ???
|| AV - Y || --> min
Решение выражается с помощью псевдообратной матрицы (Мура--Пенроуза)
A+ = pinv(A) = inv(A*A) A* (звездочкой здесь
мы обозначили транспонирование)
--------------------------------------------------
Проблема выбросов
Если в нашей выборке есть случайные ошибки, т.е. пары (x, y),
которые не подчиняются общей закономерности, то очень сильно влияют
на построенную модель.
Хотелось бы построить модель, которая минимизирует среднюю абсолютную
(линейную) ошибку:
MSE = 1/m Sum_i |f(x_i) - y_i|
(cумма модулей разности).
Но для минимизации суммы модулей разностей точного решения нет.
Минимизации методами, которые предполагают дифференцируемость
функции ошибки, препятствует недифференцируемость модуля в точке 0.
Как быть?
Можно рассмотреть компромисс: использовать функцию Хубера
{ 1/2 x^2 при |x| <= delta
H_delta(x) = {
{ delta*|x - delta/2| при |x| > delta
HuberError = 1/m Sum_i H_delta(f(x_i) - y_i)
Можно строить модель, минимизируя функцию ошибки Хубера.
Построенная модель будет значительно меньше зависеть от выбросов.
Все-таки интересно сравнить минимизацию с помощью функции Хубера
и минимизацию средней абсолютной ошибки MAE (вычисляемой при помощи модуля
разности).
Оказывается, что эту задачу можно решить методами линейного программирования.
Итак, сформулируем нашу задачу:
имеем тренировочную выборку
x1, x2, x3, ..., x_m
y1, y2, y3, ..., y_m
x_i <- R^n, y_i <- R. Надо построить линейную зависимость
y = <w, x>
При этом мы хотим минимизировать среднюю абсолютную ошибку
MAE = Sum |<w, x_i> - y_i| --> min_w
Что такое линейное программирование?
Мы имеем:
1) набор переменных, которые могут быть свободными или ограниченными
с одной или с обеих сторон
x -- свободная
y >= 0 -- ограничена снизу
z <= 3 -- ограничена сверху
1 <= t <= 12 -- ограничена с обеих сторон;
В классическом линейном программировании переменные непрерывны,
но в модуле pulp, который в Питоне используется для решения задач
линейного программирования, бывают также целочисленные переменные
и логические (бинарные);
2) набор ограничений в виде равенств или (нестрогих) неравенств:
x + y = 10
x - 2*z <= 3
5*z + 7*t >= 2
Допускаются только линейные комбинации переменных и констант;
3) одна целевая функция
2*x + 3*y - t
Мы ищем либо минимум, либо максимум целевой функции.
Нам надо переформулировать нашу задачу
MAE = Sum |<w, x_i> - y_i| --> min_w
как задачу линейного программирования.
Нам мешает модуль. Как быть???
Можно использовать простой трюк. Нам нужно избавиться от модуля ошибки
| e |
Выразим e как разность разность двух неотрицательных переменных:
e = u - v, u >= 0, v >= 0
Справедливо неравенство
|e| <= u + v
поскольку модуль суммы не превосходит суммы модулей.
Мы заменяем задачу минимизации e на минимизацию суммы u + v
|e| --> min
заменяем на
e = u - v, u >= 0, v >= 0
u + v --> min
В точке минимума либо u = 0, либо v = 0.
Если это не так и, к примеру, e > 0 и u > v > 0, то заменим u = u - v, v = 0,
в результате сумма u + v уменьшится, а равенство сохранится. Значит,
в точке минимума
|e| = u + v.
Мы свели задачу
|e| --> min
к эквивалентной задаче линейного программирования
e = u - v,
u >= 0, v >= 0
u + v --> min
В случае линейной регрессии задача переформулируется следующим образом.
Исходная задача:
Sum |<w, x_i> - y_i| --> min
(это запись в векторном виде). Через координаты:
w = (w_1, w_2, ..., w_n)
x_i = (x_i1, x_i2, ..., x_in)
y_i <- R i = 1, 2,..., m
Задача линейного программирования:
<w, x_i> - y_i = u_i - v_i
u_i >= 0, v_i >= 0, i = 1, 2,..., m
Sum_i (u_i + v_i) --> min
ПЕРЕРЫВ 10 минут до 10:45
В Питоне для решения задач линейного программирования
используется модуль pulp
from pulp import *
=============================================
29 Apr 2022
Метод опорных векторов (Support Vector Machine, SVM)
Наши объекты описываются принаками (features), которые чаще всего
представляют собой вещественные числа. Пусть у объекта n признаков,
тогда он описывается точкой n-мерного пространства R^n.
x = (x_0, x_1, ..., x_n-1) <- R^n
Пусть мы имеем 2 класса объектов. Нам нужно научиться различать
объекты этих классов. В методе опорных векторов строится линейный
классификатор:
f(x) = <w, x> - b, w <- R^n -- вектор нормали к гиперплоскости,
b -- свободный член (intercept)
b <- R
Если f(x) > 0, то мы относим объект x к первому классу,
f(x) < 0, то ко второму классу.
Если f(x) = 0, то непонятно, к какому классу отнести объект.
Уравнение
f(x) = 0
или
<w, x> - b = 0
задает гиперплоскость в n-мерном пространстве, которая разделяет
точки двух классов.
Пример: мы работаем в банке, и мы хотим написать программу,
которая определяет, сможет ли клиент вернуть кредит.
Параметры клиента описываются несколькими числами:
1) средний месячный доход;
2) общая сумма невыплаченных кредитов;
3) месячный платеж по незакрытым крелитам;
4) возраст клиента;
5) пол;
6) сколько лет является клиентом банка;
. . .
Мы строим наш классификатор по выборке: допустим, у нас имеется информация
по 1000 клиентам, для которых мы каким-то образом определили, кто является
надежным и кто нет. Мы хотим на основании этой информации построить
наш классификатор.
На тренировочных данных выполнена разметка. В методе опорных векторов
применяется СИММЕТРИЧНАЯ разметка:
Выборка:
x0, x1, x2, x3, ..., x_m-1 <- R^n -- объекты выборки
y0, y1, y2, y3, ..., y_m-1 = +-1 -- разметка выборки
y_i = { 1, если x_i принадлежит первому классу
{-1, если x_i принадлежит второму классу
Пусть мы уже построили классификатор f(x) = <w, x> - b
Тогда, если классификатор не ошибается на тренировочном объекте x_i,
то
y_i * (<w, x_i> - b) > 0
Если классификатор ошибается, то эта величина меньше 0:
y_i * (<w, x_i> - b) < 0
Как нам построить классификатор? Мы хотели бы, чтобы он не ошибался
по крайней мере на тренировочной выборке:
y_i * (<w, x_i> - b) > 0 для всех i = 0, 1, ..., m-1
Рассмотрим сначала разделимый случай.
Тогда минимальное расстояние от гиперплоскости до точек первого класса
равно минимальному расстоянию до точек второго класса, обозначим это расстояние
через h. Напомним формулу для расстояния со знаком от точки x до плоскости
<w, x> - b = 0
Расстояние
dist = (<w, x> - b)/|w|
(поскольку w -- нормаль к плоскости; единичная нормаль w/|w|)
Тогда для опорных точек выполняется точное равенство
h = +-(<w, x> - b)/|w|
Нам надо выбрать гиперплоскость так, чтобы это расстояние было максимальным:
y_i (<w, x_i> - b)/|w| >= h, i=0,1,...,m-1
h --> max_w,b
Давайте нормируем неравенства
y_i(<w, x_i> - b)/(|w|*h) >= 1
Обозначим
w' = w/(|w|*h), b' = b/(|w|*h)
получим неравенство
y_i(<w', x_i> - b') >= 1
Имеем раенство
|w'| = |w|/(|w|*h) = 1/h
То есть
h = 1/|w'|
Наша задача переписывается в виде
y_i(<w', x_i> - b') >= 1
|w'| --> min, это эквивалентно |w'|^2 --> min
Опустим штрихи, получаем задачу
y_i(<w, x_i> - b) >= 1
|w|^2 --> min_w,b
Мы получили задачу условного экстремума с ограничениями типа неравенств.
Подобные задачи решаются методами квадратичного программирования.
Мы ее сейчас решать не будем! Она была нам нужна только для того,
чтобы сформулировать нашу задачу в общем (неразделимом) случае.
Общий случай, когда точки двух классов невозможно разделить
гиперплоскостью
Пусть выполнение неравенств
y_i(<w, x_i> - b) >= 1
невозможно для всех индексов i. Тогда давайте немного ослабим эти
неравенства: введем дополнительные переменные xi_i >= 0 и разрешим
выполнение лишь ослабленных неравенств:
y_i(<w, x_i> - b) >= 1 - xi_i
Переменные xi_i называются ослабляющими переменными (slack variables),
эта как бы мера допустимой ошибки. Наша цель состоит в том,
чтобы минимизировать суммарную ошибку:
Sum xi_i --> min
Также надо минимизировать |w|^2 --> min
В методе опорных векторов используется гиперпараметр C. Итак, в общем
случае мы минимизируем целевую функцию:
|w|^2 + C*1/m * Sum xi_i --> min_w,b
y_i(<w, x_i> - b) >= 1 - xi_i
xi_i >= 0
Получили задачу квадратичного программирования: минимизировать
целевую функцию при ограничениях типа неравенств.
Оказывается, что задача условного экстремума легко сводится
к задаче безусловного экстремума. Здесь можно избавиться от
переменных xi_i. Для каждой переменной xi_i имеем два неравенства:
y_i(<w, x_i> - b) >= 1 - xi_i
xi_i >= 0
Преобразуем первое неравенство:
xi_i >= 1 - y_i(<w, x_i> - b)
xi_i >= 0
Переменная xi_i -- это мера допустимой ошибки, чем она меньше, тем лучше.
Поэтому можно положить
xi_i = max(1 - y_i(<w, x_i> - b), 0)
Подставив значение xi_i в целевую функцию, мы получаем задачу
безусловной минимизации:
|w|^2 + C*1/m * Sum_i max(1 - y_i(<w, x_i> - b), 0) --> min_w,b
Функция
H(c) = { 0, c >= 1
{ 1 - c, c < 1
называется HingeLoss
Тогда минимизируемая функция записывается в виде
|w|^2 + C*1/m * Sum_i H(y_i(<w, x_i> - b)) --> min_w,b
Мы получили задачу нахождения минимума функции.
ПЕРЕРЫВ 10 минут до 10:40
Мы реализовали на Питоне (в Jupyter-Lab) линейный метод опорных
векторов. Но часто встречается случай, когда в принципе невозможно
разделить гиперплоскостью два класса точек.
Идея: дополним признаки объектов некоторыми базовыми функциями
от признаков. Например, в простейшем случае можем добавить все
произведения (мономы) степени <= d от исходных признаков.
Например, при d = 2 вместо объекта
x = (x0, x1) <- R^2
мы рассматриваем "расширенный" объект
x' = (1, x0, x1, x0^2, x0*x1, x1^2) <- R^6
При d = 3 мы рассматриваем все мономы степени <= 3:
x' = (1, x0, x1, x0^2, x0*x1, x1^2, x0^3, x0^2*x1, x0*x1^2, x1^3) <- R^10
Если точки не разделимы гиперплоскостью в исходном пространстве,
то они с большой вероятностью будут разделимы гиперплоскостью в расширенном
пространстве (оно часто называется "спрямляющим пространством",
a отображение
phi: x --> x'
phi: R^n --> R^k, k > n
называется спрямляющим отображением).
Если мы вычислили классификатор в расширенном пространстве
f'(x') = <W, x'> - B
то классификатор в исходном пространстве будет вычисляться как
f(x) = <W, phi(x)> - B
Но разделяющая поверхность, определяемая уравнением
<W, phi(x)> - B = 0
уже не будет плоскостью. Таким образом, мы построим
нелинейный классификатор.
=========================================================
6 May 2022
Ядровый метод опорных векторов
Классический линейный метод опорных векторов
Даны объекты двух классов
x1, x2, ..., x_m -- тренировочная выборка
x_i <- R^n
Разметка симметричная
y1, y2, ..., y_m, y_i = +-1
y_i = { 1, если x_i принадлежит первому классу
{ -1, если x_i принадлежит второму классу
Требуется построить линейный классификатор:
f(x) = <w, x> - b, w <- R^n, b -- свободный член, b <- R
Если f(x) > 0, то мы относим x к первому классу,
f(x) < 0, то ко второму.
Если f(x) = 0, то ситуация неопределенная.
Уравнение
f(x) = 0
задает нам гиперплоскость коразмерности 1 в R^n
<w, x> - b = 0
Эта гиперплоскость разделяет пространство на 2 полупространства,
мы хотим, чтобы точки первого класса лежали бы в верхнем полупространстве
<w, x> - b > 0,
а точки второго класса в нижнем полупространстве
<w, x> - b < 0.
Заметим, что точка x_i лежит в правильном полупространстве <=>
y_i (<w, x_i> - b) >= 0
Расстояние cо знаком от точки x до плоскости
<w, x> - b = 0
равно
(<w, x> - b)/||w||
Мы получаем, что для опорных точек расстояние h
до разделяющей гиперплоскости равно:
y_i (<w, x_i> - b) / ||w|| = h
Для остальных точек (в разделимом случае)
y_i (<w, x_i> - b) / ||w|| > h
Мы получаем задачу:
y_i (<w, x_i> - b) / ||w|| >= h, i=1,2,...,m
h --> max_w,b
Давайте избавимся от h
y_i (<w, x_i> - b) / (|w|*h) >= 1
y_i (<w/(|w|*h), x_i> - b/(|w|*h)) >= 1
Обозначим
w/(|w|*h) = w', b/(|w|*h) = b'
Получим
y_i (<w', x_i> - b') >= 1
h --> max
Выразим h через |w'|
w' = w/(|w|*h)
|w'| = |w|/(|w|*h) = 1/h
h = 1/|w'|
Таким образом,
h --> max
эквивалентно
|w'| --> min
или, что эквивалентно,
|w'|^2 --> min
Опуская штрихи, получаем задачу:
y_i (<w, x_i> - b) >= 1, i=1, 2, ..., m
|w|^2 --> min_w,b
Разделяющая гиперплоскость находится как решение этой оптимизационной задачи.
Она решается методами квадратичного программирования.
В неразделимом случае эта задача неразрешима.
НЕРАЗДЕЛИМЫЙ СЛУЧАЙ
Все неравенства
y_i (<w, x_i> - b) >= 1, i=1, 2, ..., m
выполнить невозможно. Давайте ослабим их, введя ослабляющие переменные
ksi_i >= 0 для каждого их неравеств:
y_i (<w, x_i> - b) >= 1 - ksi_i, i=1, 2, ..., m
Переменные ksi_i называются slack variables.
Значение ksi_i -- это как бы штраф за неточное выполнение неравенства.
Наша цель -- минимизировать суммарный штраф + максимизировать ширину полосы,
разделяющей два класса точек. Для балансирования этих двух целей
используется гиперпараметр C > 0. Получаем задачу
условной оптимизации:
|w|^2 + C*(1/m) * Sum_i ksi_i --> min_w, b
y_i (<w, x_i> - b) >= 1 - ksi_i,
ksi_i >= 0, i=1,2,...,m
На прошлых лекциях мы избавились от переменных ksi_i и свели задачу к
безусловной оптимизации:
y_i (<w, x_i> - b) >= 1 - ksi_i,
ksi_i >= 0
Перенесем ksi_i в левую часть:
ksi_i >= 1 - y_i (<w, x_i> - b),
ksi_i >= 0
Поскольку наша цель -- сделать ksi_i как можно меньше, то можно взять
ksi_i = max(1 - y_i (<w, x_i> - b), 0)
Мы получаем задачу безусловной минимизации:
|w|^2 + C*(1/m) * Sum_i max(1 - y_i (<w, x_i> - b), 0) --> min_w, b
Функция
H(c) = max(1 - c, 0)
называется Hinge Loss. Задача переписывается в виде:
|w|^2 + C*(1/m) * Sum_i H(y_i (<w, x_i> - b)) --> min_w, b
Для минимизации функции можно использовать модуль Питона
from scipy.optimize import minimize
res = minimize(f, x0)
res.x -- точка минимума
Прошлый раз был рассмотрен нелинейный метод разделения.
Идея: дополним признаки объекта дополнительными признаками,
которые вычисляются как некоторые базовые функции от исходных
признаков.
x = (x1, x2, ..., x_n) <- R^n
Рассмотрим набор фунций
phi_1(x1, x2, ..., x_n), phi_2(x1, x2, ..., x_n), ...,
phi_k(x1, x2, ..., x_n)
Обычно k > n
Заменим точку x на x'
x' = (phi_1(x), phi_2(x), ..., phi_k(x)) <- R^k
phi: x --> x'
phi: R^n --> R^k -- спрямляющее отображение.
Можно надеяться, что в спрямляющем пространстве R^k наши точки x'_i
уже будут разделяться гиперплоскостью. Мы находим линейный классификатор
в расширенном пространстве
F(x') = <W, x'> - B
Тогда в исходном пространстве классификатор выражается формулой
f(x) = <W, phi(x)> - B
Получаем нелинейный классификатор.
В простейшем случае в качестве спрямляющего отображения phi
можно рассмотреть все мономы степени <= d от исходных признаков.
Например, пусть n = 2, x = (x0, x1), d = 2,
phi(x) = (1, x0, x1, x0^2, x0*x1, x1^2) <- R^6
При d = 3 имеем
phi(x) = (1,
x0, x1,
x0^2, x0*x1, x1^2,
x0^3, x0^2*x1, x0*x1^2, x1^3) <- R^10
Проблема: мы сводим исходную задачу в пространстве R^n к пространству R^k,
k > n, в результате получаем существенное замедление вычислений.
Оказывается, что можно переходить в расширенное (спрямляющее) пространство
и при этом скорость вычислений существенно не уменьшается.
Для этого используется ядровый трюк.
Определение (из функционального анализа).
Функция s = K(x, y), x, y <- R^n, s <-R -- число
называется ядром, если существует отображение
phi: R^n --> H -- Гильбертово пространство
такое, что
K(x, y) = <phi(x), phi(y)>
Ядро -- это скалярное произведение образов векторов x, y в
некотором расширенном пространстве (оно может быть даже бесконечномерным).
ВАЖНО!!!
Для вычисления значения ядра K(x, y) нам не нужно знать отображение
phi и реально вычислять образы векторов phi(x), phi(y), поскольку
ядро обычно задано просто какой-то формулой.
Например, Гауссовское ядро RBF (Radial Basic Function)
K(x, y) = exp(-||x - y||^2/2*sigma^2)
ПЕРЕРЫВ 10 минут до 10:50 (Moscow)
Вернемся к нашей задаче. Цель -- переформулировать ее так,
чтобы в ней использовались только скалярные произведения
исходный векторов <x_i, x_j>. Получим задачу в исходном пространстве.
Если мы теперь заменим скалярные произведения на значение ядра
<x_i, x_j> <---> K(x_i, x_j)
то это будет эквивалентно переходу к задаче в расширенном пространстве
для векторов phi(x_i), i=1,2,...,m. Но при этом вычисления не будут
усложняться (просто вместо скалярного произведения всюду будем
использовать значение ядра). В результате мы получим классификатор:
f(x) = K(w, x) - b
Итак, вернемся к задаче условной оптимизации:
целевая функция
|w|^2 + C*(1/m) * Sum_i ksi_i --> min_w, b
ограничения типа неравенств
y_i (<w, x_i> - b) >= 1 - ksi_i,
ksi_i >= 0, i=1,2,...,m
Для задачи условного экстремума имеется теорема Карруша-Куна-Такера,
которая обобщает теорему Лагранжа об условном экстремуме.
Теорема Лагранжа:
f(x) --> min (или мах) локальный
x <- R^n
h1(x) = 0
h2(x) = 0
...
h_k(x) = 0
В точке экстремума x* выполняется условие равенства нулю
градиента лагранжиана:
L(x, lambda_1, lambda_2, ..., lambda_k) =
f(x) + lambda_1 h1(x) + ... + lambda_k h_k(x)
d/dx (L(x*)) = 0 градиент лагранжиана по x равен нулю
для некоторых коэффициентов lambda.
В теореме Лагранжа все ограничения типа равенств.
Теорема Куна-Такера обобщает теорему Лагранжа, когда
есть ограничения типа (нестрогих) неравенств.