Факультет фундаментальной
физико-химической инженерии
2021-22 учебный год
Химики, 2-й курс, весенний семестр
12 Feb 2022
Язык программирования Python3
-----------------------------
Почему Питон?
1. Он очень легко учится. Из-за этого он стал очень популярен в мире.
2. В Питоне нет компилятора, язык интерпретируемый, можно с Питоном
работать в интерпретаторе и использовать Питон в качестве
очень удобного и мощного калькулятора.
3. За счет простоты и популярности Питон приобрел массу полконников
(создатель языка Гвидо Россум), из называют питонистами.
В результате написана огромное количество программ на Питоне
и множество модулей для Питона. Благодаря этому разрабатывать
свои приложения на Питоне довольно просто. Например, почти весь
код, связанный с искусственным интеллектом, написан на Питоне
(нейросети, машинное обучение и т.п.) -- вернее, Питон использует
модули, написанные на C++, и собсвенно на Питоне пишется только
самый верхний уровень.
Питон -- это объектно-ориентированный язык с динамической типизацией.
В Питоне вообще нет описаний переменных!
C++:
std::vector a(5) = { 1.5, 2, 3.7, 4, 6 };
Python:
a = [ 1.5, 2., 3.7, 4., 6. ]
Как установить Питон?
MS Windows -- установить Python 3.7
+ можно установить Jupyter-Lab
Linux -- просто устанавливаете python3 c помощью apt
sudo apt install python7
Первая программа на Питоне: Hello, World!
print("Hello, Python!")
Совсем примитивно.
Давайте напишем программу, которая печатает квадраты первых 20
натуральных чисел.
Пишем в ДУРНОМ СТИЛЕ (не вздумайте так писать!!!)
n = 1
while n <= 20:
print(n, n*n)
n += 1
ВНИМАНИЕ!!!
Дурной стиль -- иметь в программе код на глобальном уровне
(вне функций), который выполняется безусловно.
Это препятствует использованию файла с программой как модуля
в других программах Питона. Код на глобальном уровне должен выполняться
только тогда, когда программа выполняется как скрипт.
В этом случае специальная переменная "Имя модуля"
__name__
имеет значение "__main__". Иначе она равна основному имени
файла с Питон-программой
(имена файлов на Питоне должны состоять из латинских букв, цифр
и символа подчеркивания и не могут начинаться с цифры!!!).
def printSquares(upperBound = 20):
n = 1
while n <= upperBound:
print(n, n*n)
n += 1
if __name__ == "__main__":
printSquares()
Даже лучше весь код, выполняемый в случае запуска программы
как скрипт, поместить в функцию main()
def printSquares(upperBound = 20):
n = 1
while n <= upperBound:
print(n, n*n)
n += 1
def main():
print("Squares of the first n integers")
while True:
try:
n = int(input("n: "))
break
except ValueError:
print("Incorrect input")
printSquares()
if __name__ == "__main__":
main()
Мораль: весь код должен быть внутри функций.
Функция, которая выполняется, когда программа запускается
как скрипт, должна называться main()
или как-нибудь еще. Две последние строки в файле -- это всегда
if __name__ == "__main__":
main()
Пусть мне нужно не напечатать что-то, а получить список квадратов.
Список в Питоне -- это на самом деле динамический массив,
аналог std::vector в C++. Элементы списка указываются в
квадратных скобках:
smallPrimes = [2, 3, 5, 7, 11, 13, 17, 19]
Элементы списка обозначаются как smallPrimes[i],
где индекс i -- это целое число от 0 до len(smallPrimes) - 1
Пусть a -- это список. Первый элемент списка -- это a[0],
последний элемент -- это a[len(a) - 1]
Но Гвидо Россум сделал обращение к последним элементам списка
очень удобным: последний элемент -- это a[-1]
a[-2] -- это предпоследний элемент и т.д.
Помимо списков, в Питоне есть tuple (кортеж), который
обобщает понятие пары элементов, тройки элементов, четверки и т.п.
pair пара
triple тройка
quadruple четверка
. . .
tuple n-ка энка кортеж
Чем отличается tuple от списка?
Элементы списка можно менять, элементы tuple менять уже нельзя
(immutable objects).
Первая тема -- это примения теории чисел в криптографии.
Имейте в виду, что Питоне3 операция деления целых чисел
возвращает вещественный результат! (Происходит потеря информации.)
Операция целочисленного деления обозначается //
>>> 5//2
2
>>> 5/2
2.5
Комментарий в Python обозначается #,
или многострочный комментарий
'''
Multiline
comment
'''
Способ задания списка, характерный для Питона,
который питонисты применяют на каждом шагу:
List Comprehension ("постижение списка" -- бессмыслица???)
"Умное задание списка" или "Математическое задание списка".
Как в математике задать список всех нечетных чисел, меньших 100?
S = { n: n <- 1, 2, ..., 100, n нечетно }
В Питоне:
s = [ n for n in range(1, 101) if n%2 != 0 ]
В List Comprehension мы сначала указываем выражение,
зависящее от параметров, затем указываем диапазон
изменения каждого параметра, а также можем задать фильтр
(предикат), который отбирает элементы, включаемые в список.
Пример: квадраты натуральных чисел
squares = [(n, n*n) for n in range(1, 21)]
Мы написали функцию isPrime(m), которая возвращает True для
простых чисел m. Как получить список всех простых чисел,
меньших 100? Используем list comprehension:
smallPrimes = [ m for m in range(100) if isPrime(m) ]
Пример 2: найти все простые числа-близнецы, меньшие 1000:
primeTwins = [(m, m+2) for m in range(999)
if isPrime(m) and isPrime(m + 2)]
-----------------------------------------------------------
19 Feb 2022
Прошлый раз была рассказана конструкция List Comprehension
(умное задание списка)
[ выражение от параметров for диапазон изменения параметров
if фильтр ]
[ m for m in range(1, 101) if m%3 == 0 and m%7 == 0 ]
Задача: найти все пифагоровы тройки (a, b, c) в диапазоне до 100
a^2 + b^2 = c^2
числа a, b, c целые, 1 <= a < b < c <= 100,
gcd(a, b) = 1
(если тройка (a, b, c) пифагорова, то тройка (a*m, b*m, c*m)
тоже пифагорова
Например, (3, 4, 5) пифагорова тройка (египетский треугольник)
то (6, 8, 10),
(9, 12, 15),
. . .
[ (a, b, c) for a in range(1, 101)
for b in range(a+1, 101)
for c in range(b+1, 101)
if a*a + b*b == c*c ]
Задача
Найти все счастливые билеты длины n, где n > 0 -- четное число
(сумма цифр первой половины номера билета равна сумме цифр
торой половины).
Будем представлять номер билета как список цифр от 0 до 9
Домашнее задание 1
Написать функцию для числа счастливых билетов длины n,
которая будет работать разумное время для n <= 12.
Домашнее задание 2
Написать функцию, проверяющую, является ли число m полусовершенным.
Число m совершенное, если оно равно сумме своих собственных
делителей.
6 = 1 + 2 + 3
28 = 1 + 2 + 4 + 7 + 14
Число m называется ПОЛУсовершенным, если оно равно сумме
какого-то ПОДМНОЖЕСТВА своих собственных делителей.
Как эту задачу решать?
1) находим список d всех собственных делителей числа m;
2) решаем задачу о рюкзаке napsack(d, m).
Задача о рюкзаке.
Имеется рюкзак объема m и список вещей, в содержащий объемы
каждой вещи. Можно ли какими-то вещами из этого списка
(не обязательно всеми) заполнить рюкзак доверху?
[1, 2, 5, 3, 7, 11] m = 10
Можно ли из этого списка выбрать числа,
сумма которых равна 10?
2, 5, 3
3, 7
1, 2, 7
Эта задача NP-полная и для нее не существует хорошего
решения, кроме как полного перебора всех подмножеств.
Домашнее задание 3.
Дана перестановка s чисел 1, 2, ..., n,
представленная списком длины n.
Написать функцию nextPermutation(s),
которая в том же самом списке получает следующую
перестановку в лексикографическом порядке.
Функция возвращает True, если это можно сделать,
и False, если нельзя.
Идея
1 2 3 4
1 2 4 3
1 3 2 4
1 3 4 2
1 4 2 3
1 4 3 2
2 1 3 4
2 1 4 3
2 3 1 4
2 3 4 1
2 4 1 3
. . .
4 3 2 1
=========================================
26 Feb 2022
Введение в Python3
Элементы теории чисел и кодирование с открытым ключом
В Питоне поддерживаются целые числа неограниченного размера.
Элементы теории чисел
Главный объект -- это кольцо Z_m вычетов по модулю m.
Целые числа Z -- кольцо (есть операции +, *)
Назовем два числа x1, x2 эквивалентными
x1 ~ x2
если m | (x1 - x2) или если (x1 - x2) <- mZ
Это обычно обозначается так:
x1 == x1 (mod m)
читается как x1 равно x2 по модулю m
Z/mZ -- фактор-кольцо Z_m
Кольцо целых чисел бесконечно, а кольцо Z_m конечное,
в нем m элементов.
Компьютер всегда работает с элементами кольца вычетов Z_m,
а не с целыми числами (которых бесконечно много).
Как представляются элементы Z_m?
Неотрицательная система остатков
Z_5 = { 0, 1, 2, 3, 4 }
_ _ _ _
2*4 == 8 == 3 (mod 5)
Симметричная система остатков:
Z_5 = { -2, -1, 0, 1, 2 }
Z_8 = { -4, -3, -2, -1, 0, 1, 2, 3 }
-4 == 4 (mod 8) 4 + 4 == 0 (mod 8)
Целые числа в компьютере (в языке C) -- это элементы Z_m,
где m = 2^32 x = 2^31 == -2^31 (mod m)
x + x == 0 (mod m)
x + y = ???
(x + y)%m
При выполнении операций в кольце Z_m надо всегда от результата
оставлять только остаток от деления на m. 3, 8, -2
Число -1 и число m-1 -- это одно и то же в кольце Z_m
В частности, в Z_5
-1 == 4 (mod 5)
Самая знаменитая (и самая полезная) теорема в теории чисел -- это
Малая теорема Ферма.
Пусть p -- простое число. Тогда для всякого числа b != 0 (mod p)
выполняется
b^(p-1) == 1 (mod p)
(любой ненулевой элемент кольца Z_p является корнем p-1 -ой степени из
единицы).
Иногда используют такую формулировку: для всякого b
b^p == b (mod p)
Возведение в степень p -- это тождественное отображение кольца Z_m
b^p = b * b^(p-1) == b
B даже такую:
для всякого s >= 0 и для любого b
b^(s*(p-1) + 1) == b (mod p)
Большая теореме Ферма:
существуют ли n > 2, целые a, b, c такие, что
a^n + b^n = c^n ?
В отличие от Малой теоремы Ферма, Большая на практике не используется.
Теорема Эйлера
Рассмотрим группу обратимых элементов по умножению
кольца Z_m, где m -- не обязательно простое число.
U_m лежит в Z_m
Элемент x обратим в кольце Z_m <==> gcd(x, m) = 1
x взаимно прост с m
x обратим, если найдется y: x*y == 1 (mod m)
Основная теорема арифметики:
рассмотрим два числа m, n, не равных одновременно не равных нулю
(m != 0 или n != 0)
Тогда d = gcd(m, n) выражается в виде
d = u*m + v*n
Число d вычисляется с помощью расширенного алгоритма Евклида.
Просто gcd(m, n) вычисляется следующим образом:
(m, n) --> (n, r) --> . . . -> (d, 0)
r -- остаток от деления m на n
m = q*n + r, |r| < |n|
r = m - q*n
gcd(m, 0) = m
def gcd(m, n):
while n != 0:
(m, n) = (n, m%n)
return m
Построим расширенный алгоритм Евклида
Дано: m, n (хотя бы одно из этих чисел отлично от нуля)
Надо: вычислить d = gcd(m, n) и числа u, v такие, что
d = u*m + v*n
Идея: применим схему построения цикла с помощью инварианта.
Будем применять обычный алгоритм Евклида, параллельно
(a, b)
a = m, b = n
(a, b) --> (b, r) --> ... --> (x, 0)
При этом на каждом шаге мы храним коэффициенты выражения
чисел a, b через исходные числа m, n
Инвариант цикла:
gcd(a, b) = gcd(m, n)
a = u1*m + v1*n
b = u2*m + v2*n
Начальные значения:
a = m, b = n
u1 = 1, v1 = 0
u2 = 0, v2 = 1
Действие в цикле
Делим с остатком: находим q, r такие, что
a = q*b + r
(a, b) --> (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)
В конце b = 0, d = a (наибольший общий делитель),
u = u1, v = v1
Зачем нужен Расширенный алгоритм Евклида?
Он нужен для нахождения обратного элемента к x в кольце вычетов Z_m
Дано: x, m
Надо: найти y такой, что x*y == 1 (mod m)
Решение: применив Расширенный алгоритм Евклида к паре (x, m),
найдем (d, u, v)
d = u*x + v*m
Если d != 1, то x необратим
Если d = 1, то y = u -- обратный элемент
1 = u*x + v*m == u*x (mod m), так как v*m == 0
Почему может быть нужен обратный элемент?
Кодировка
фиксируем большое число m и е, взаимно простое с m
Рассмотрим отображение (кодировка)
E: x --> x*e
===========================================================
7 Mar 2022
Было рассказано: что такое кольцо вычетов по модулю m
Малая теорема Ферма:
Пусть p -- простое число, и пусть b -- любое число, которое не
делится на p:
b != 0 (mod p)
Тогда b^(p-1) == 1 (mod p).
Варианты:
Пусть p -- простое, b -- любое число. Тогда
b^p == b (mod p)
(возведение в степень p -- это тождественное отображение в кольце
вычетов по модулю p).
И еще один вариант:
Пусть p -- простое, b -- любое число, k >= 0 -- тоже любоу целое неотрицательное
число. Тогда
b^(k*(p - 1) + 1) == b (mod p)
(возведение в степень k*(p - 1) + 1 -- это тождественное отображение в кольце
вычетов по модулю p).
Доказательство.
При b == 0 очевидно.
При b != 0 (mod p)
b^(k*(p - 1) + 1) = b^(k*(p - 1)) * b =
= (b^(p - 1))^k * b == (1)^k * b == b (mod p)
Китайская теорема об остатках
Пусть m = m1 * m2 * ... *m_k, gcd(m_i, m_j) = 1
Примеры: m = 105 = 3*5*7
m = 100 = 4*25
m = p1^e1 * p2^e2 * ... * pk^ek
Рассмотрим Z_105 105 = 3*5*7
56 = (2, 1, 0)
19 = (1, 4, 5)
+
75 = (0, 0, 5)
17 = (1, 2, 3)
10, 17, 24, 31, 38, 45, 52, 59, 66, 73, 80
10 = (1, 0, 3)
17 = (2, 2, 3)
52 = (1, 2, 3) Загайнова +
Китайская теорема об остатках позволяет свести изучение кольца Z_m
для произвольного m к изучению кольца Z_p^e, где p простое.
Теорема Эйлера
Обозначим через phi(m) -- число обратимых элементов в кольце Z_m.
Пример: чему равно phi(30)?
~
Z_30 = Z_2 + Z_3 + Z_5
(1, 1, 1)
(1, 1, 2)
(1, 1, 3)
(1, 1, 4)
(1, 2, 1)
(1, 2, 2)
(1, 2, 3)
(1, 2, 4)
phi(30) = phi(2*3*5) = (2-1)*(3-1)*(5-1) = 8
Для всякого b != 0 и для всякого m = p1^e1 * p1^e2 * ... * p_k^e_k
phi(m) = (p1-1)*p1^(e1-1) * (p2-1)*p2^(e2-1) * ... * (p_k-1)*p_k^(e_k-1)
В частности, для числа m, свободного от квадратов (т.е. не делящегося на
квадраты любых числел)
m = p1 * p2 * ... * p_k
функция Эйлера
phi(m) = (p1 - 1)*(p2 - 1)*...*(p_k - 1)
Пример:
phi(100) = phi(2^2 * 5^2) = (2-1)*2^1 * (5-1)*5^1 = 2*20 = 40
Теорема Эйлера
Для всякого b, взаимно простого с m, gcd(b, m) = 1
выполняется
b^phi(m) == 1 (mod m)
Форма теоремы Эйлера:
пусть m свободно от крадратов, m = p1*p2*p3*...*p_k, где p_i -- простые числа.
Пусть s >= 0. Тогда для всякого b
b^(s*phi(m) + 1) == b (mod m)
Схема кодирования с открытым ключом RSA
(MIT professors Rumley, Shamir, Adleman).
Пусть p, q -- два больших простых числа, например, 200-значных десятичных.
Пусть m = p*q
Число m открыто, разложение m на множители секретно (числа p и q держатся в
секрете).
Вычислим phi(m) = (p - 1)(q - 1)
Возьмем любой элемент e, обратимый в кольце вычетов Z_phi(m)
Вычисляем элемент d, обратный к e в кольце Z_phi(m), используя
расширенный алгоритм Евклида.
d*e == 1 (mod phi(m))
Открытый ключ: пара (m, e)
Секретный ключ: пара (m, d)
Отображение шифрования:
E: Z_m --> Z_m (Encoding)
E: x --> x^e (mod m)
Обратное отображение (расшифровка):
D: Z_m --> Z_m (Decoding)
D: y --> y^d (mod m)
Докажем, что D*E = E*D = 1 -- тождественное отображение.
x -- исходное сообщение
y == x^e (mod m) -- зашифрованное сообщение
z == y^d (mod m) = (x^e)^d = x^(e*d) = x^(1 + s*phi(m)) == x (mod m)
по теореме Эйлера
e*d == 1 (mod phi(m))
e*d = 1 + s*phi(m)
Пример:
p = 997 q = 1019
m = p*q = 1015943
phi = (p - 1)*(q - 1) = 1013928
e = 103
d = invmod(e, phi) = 255943
Открытый ключ: пара (m = 1015943, e = 103)
Серкетный ключ: пара (m = 1015943, d = 255943)
Исходное сообщение:
x = 12345
Зашифрованное сообщение:
y = x^e (mod m) = 783688
Расшифруем его:
z = y^d (mod m) = 12345
Классические шифры были симметричными:
зная ключ шифрования, млжно было легко расшифровать сообщение.
Кодирование с открытым ключом асимметричное (public key encryption).
Прямая схема: канал от всех к одному.
Еще более важна обратная схема: канал от одного ко всем.
Автор шифрует сообщение с помощью своего обратного ключа:
y = D(x)
Расшифровать его может любой (поскольку все знают открытую процедуру E)
E(y) = E(D(x)) = x
Но создать y может только обладатель секретной процедуры.
Таким образом, он удостьоверяет свою личность.
На том же принципе может быть основана электронная подпись.
Мы посылаем открытый документ и к нему посылаем дополнительно
электронную подпись, удостоверяющую его подлинность.
1. Считаем длину документа len.
2. Вычисляем hash-функцию от содержимого документа.
3. Добавляем номер транзакции и текущее время.
4. Ко всему этому применяем свою секретную процедуру D.
Получаем подпись документа.
Для проверки получатель документа сначала вычисляет его длину
и hash-функцию от его содержимого.
Затем он расшифровывает эл. подпись, используя всем известную
открытую процедуру E и сравнивает длину и hash-функцию, полученную
из эл. подписи, с реальной длиной и hash-функцией документа.
==========================================================
12 Mar 2022
Применение теории чисел в программировании.
Кодирование с открытым ключом
m = p*q
e -- обратимо в кольце Z_phi(m)
phi(m) = (p-1)(q-1)
d -- обратный элемент к e
d*e == 1 (mod phi(m))
(m, e) -- открытый ключ
(m, d) -- секретный ключ
Шифрование: E: Z_m --> Z_m
x --> x^e (mod m)
Дешифровка: D: Z_m --> Z_m
y --> y^d (mod m)
Это две взаимно обратные процедуры:
DE = 1
ED = 1
------------------------------------
Задача, возникающие в связи со схемой RSA
1. Для генерации пары ключей надо уметь генерировать
большие простые числа (примерно 200-значные десятичные).
Это мы можем делать, а разложить на множители 400-значное число
невозможно современными алгоритмами.
2. Для взлома ключа надо уметь раскладывать число на множители
(факторизация). Чтобы быть уверенным в том, что код невозможно
взломать, надо изучить возможные алгоритмы факторизации.
Тест простоты
Дано: число m
Надо: определить, является ли m простым числом.
Теорема Чебышева
Количество простых чисел <= N ~ N/ln(N)
В отрезке длины ln(N) в среднем содержится одно простое число.
В качестве теста простоты можно попробовать использовать
Малую теорему Ферма: m простое => для всякого b != 0 (mod m)
b^(m-1) == 1 (mod m)
В частности,
2^(m-1) == 1 (mod m)
Это необходимое условие простоты. Древние греки ошибочно считали,
что это и достаточное условие:
2^(m-1) == 1 (mod m) ==> m простое.
ЭТО НЕВЕРНО! Контрпример привел французский математик Сарруз в
XIX веке:
m = 341 = 31*11
2^340 == 1 (mod 341)
2^340 = (2^10)^34 = 1024^34 == 1^34 == 1 (mod 341)
1024 = 3*341 + 1 == 1 (mod 341)
Это называется тест Ферма для основания 2.
Но 341 не выдерживает тест Ферма для основания 3:
3^340 == 56 (mod 341)
Оказывается, что есть такие парадоксальные числа m, которые ведут себя
как простые в Малой теореме Ферма для всех оснований b, взаимно простых
c m:
для любого b: gcd(b, m) = 1
выполняется утверждение Малой теоремы Ферма:
b^(m-1) == 1 (mod m)
Первое такое число 561 нашел американский математик Р.Каймайкл в 1910 г.
В его честь такие числа называются кармайкловыми. Их бесконечно много:
561, 1105, 1729, ...
Для кармайкловых чисел тест простоты Ферма не работает.
Однако можно его чуть модифицировать, чтобы он работал для любых чисел.
Этот тест был найден несколькими математиками в 60-x годах прошлого века.
Тест Миллера--Рабина
Это вероятностный тест. Он исползует датчик случайных чисел для
генерации основания b. Тест может выдавать 2 ответа:
1) число m составное (и тогда выдается доказательство этого факта);
2) не знаю.
Причем для составного числа m вероятность ответа "не знаю" <= 1/4.
Если при 20 запусках мы получили 20 ответов "не знаю", то, если число m
было бы составным, то вероятность 20-ти таких ответов <= (1/4)^20.
То есть, скорее всего, число m простое, но мы не получаем доказательства
этого факта. Однако для практических целей этого достаточно (вероятность
ошибки можно сделать сколь угодно маленькой).
Как работает тест Миллера-Рабина
В тесте Ферма мы берем случайное основание b и возводим его в степень
m-1 в кольце Z_m:
b^(m-1) (mod m)
Поскольку m нечетное число, то m-1 -- четное число, вытащим из него
максималбную степень двойки:
m = t*2^s, где t нечетное.
Например, пусть m = 561, m-1 = 560 = 35 * 2^4.
Для вычисления степени m-1 сначала возведем b в степень t,
а затем s раз возведем полученное число в квадрат (все вычисления
в кольце Z_m).
x0 == b^t (mod m),
x1 == x0^2 (mod m),
x2 == x1^2 (mod m),
...
x_s == x_s-1^2 == b^(m-1) (mod m)
Итак, получили последовательность
x0, x1, x2, ..., x_s
Утверждение
если последовательность имеет вид либо
*, *, *, ..., *
либо
*, *, *, ..., *, 1, 1, ..., 1
где * != +-1 (mod m),
то число m составное.
Ответ "не знаю" выдается в случаях, когда последовательность
имеет вид либо
1, 1, 1, ..., 1
либо
*, *, ..., *, -1, 1, ..., 1
======================================================
Предложение:
тест не врет, когда он говорит, что m составное.
если последовательность имеет вид
*, *, *, ..., *
то не выполняется Малая теорема Ферма.
Если последовательность имеет вид
*, *, *, ..., *, 1, 1, ..., 1
====
то имеем x_i != +-1, x_i^2 == 1 (mod m)
Если бы m было простым, то Z_m было бы полем, а в поле только 2
квадратных корня из 1: 1 и -1.
(Теорема Безу: число корней многочлена не превышает его степени.
Корень квадратный из 1 -- это корень многочлена
f(x) = x^2 - 1
).
Второе утверждение, что вероятность последовательностей вида
либо
1, 1, 1, ..., 1
либо
*, *, ..., *, -1, 1, ..., 1
для состваного m не превышает 1/4, тоже легко доказывается
(идея: что любая подгруппа группы имеет не больше n/2 элементов).
Пример:
применим тест Миллера-Рабина для число 561
Берем случайное основание b = 2
560 = 35 * 2^4
x0 = 2^35 == 263 (mod 561),
263^2 == 166,
166^2 == 67,
67^2 == 1
263, 166, 67, 1, 1
=> число m составное
===================================================
Алгоритмы факторизации
Есть два очень простых и очень красивых алгоритма факторизации Полларда,
они часто называются методами Монте-Карло.
Интересно, что программа получается совсем простая, зато объяснить,
почему она работает, довольно сложно (т.е. алгоритм нетривиальный
и требует некоторой математической культуры).
rho-алгоритм Полларда основан на поиске цикла в
рекуррентной последовательности.
Дано: число m
Надо: найти нетривиальный делитель числа m.
Работаем в кольце вычетов Z_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+1 = f(x_i), ...
x0, x1, x2, x3, x4, ... -- орбита элемента x0
Так как Z_m конечно, то последоватленость зациклится.
x_i == x_i+d (mod m)
Пусть p | m -- нетривиальный делитель числа m. Мы его не знаем,
однако теоретически можно рассмотреть те же элементы по модулю p:
_ _ _ _ _ _
x0, x1, x2, x3, x4, ... x_i == x_i (mod p) -- элементы Z_p
Поскольку Z_p меньше, чем Z_m, то, скорее всего, в последовательности
_
{ x_i } зацикливание произойдет раньше, чем в { x_i }
x_i == x_j (mod p), j > i
x_i != x_j (mod m)
то есть
p | (x_i - x_j)
m не делит (x_i - x_j)
==> p | gcd(x_i - x_j, m), m не делит gcd(x_i - x_j, m) ==>
d = gcd(x_i - x_j, m) -- нетривиальный делитель числа m.
Как искать цикл в последовательности?
x0 <---> x1
x1 <---> x3
x2 <---> x5
x3 <---> x7
x4 <---> x9
. . .
x_i <--> x_2i+1
Оценка времени работы такого алгоритма получается как средняя длина
цикла в орбите произвольного элемента при полиномиальном
отображении f: Z_p --> Z_p Это p^(1/2).
В свою очередь, p^2 <= m ==>
t = O(m^(1/4))
Например, для 20-значного десятичного числа количество шагов порядка 10^5 =
10,000.
=======================================
19 Mar 2022
p-1 алгоритм факторизации Полларда
Дано: целое число m
Надо: найти нетривиальный делитель d | m, d != 1, d != m
Идея: воспользуемся Малой теоремой Ферма
Пусть p -- простое число; пусть b != 0 (mod p).
Тогда b^(p-1) == 1 (mod p)
Пусть m -- составное число, и пусть p | m -- простой делитель.
Предположим, что число p-1 -- гладкое (smooth).
В определении гладкости участвует константа N (например, N = 100_000_000).
Число k гладкое, если
k = q1^e1 * q2^e2 * ... * q_l^e_l, причем q_i^e_i <= N
Так как число p-1 -- гладкое, то N! делится на p-1
p-1 = m1 m2 m3 ... *m_l, m_i взаимно простые и m_i <= N
N! = 1*2*3*4*...*N -- все m_i содержатся в этом произведении
Возьмем случайное число b: gcd(b, m) = 1 (если не 1, то это уже
нетривиальный делитель).
Тогда b^(N!) == 1 (mod p), поскольку N! = (p-1)*s
b^(N!) = (b^(p-1))^s == (1)^s == 1 (mod p)
Получили, что p | b^(N!) - 1 ==> p | gcd(b^(N!) - 1, m)
В алгоритме мы берем число b и последовательно его возводим по модулю p
в степень 2, то, что получилось, в степень 3, потом в степень 4, ...
до N. На каждом шаге вычисляем gcd(b - 1, m). Если повезет, то на каком-то
шаге этот НОД будет нетривиальным.
Давайте напишем фунцию на Питоне.
. . .
Имеется вариант этого алгоритма за авторством Ленстра, который работает
на эллиптических кривых.
y^2 = x^3 + a*x + b
На точках элл. кривой можно определить операцию +,
относительно этой операции точки элл. кривой превращаются в
абелеву группу. Важно, что координаты суммы двух точек p+q
вычисляются с помощью четырех арифметических операций +, -, *, /.
Следовательно, можно определить элл. кривую над любым полем,
в том числе над конечным полем Z_p.
Малая теорема Ферма: порядок элемента группы ненулевых элементов
поля Z_p делит порядок группы, т.е. p - 1
x^e = 1, то минимальное такое e называется порядком.
x^e = x*x*x*...*x
e раз
Когда операция +, то аналогом возведения в степень является умножение
на целое число:
x*e = x+x+x+...+x
e раз
Для группы элл. кривой порядок элемента делит порядок группы.
Если k -- порядок группы, то
b*k = 0 в группе элл. кривой, т.е. ее бесконечной точке
b = (x, y, z)
b*k = ( , *, 0) (mod p) так как кривая рассматривается над Z_p
Алгоритм:
берем случайную точку на элл. кривой и последовательно
ее умножаем на 2, 3, 4, 5, 6, ..., N
действия производим в Z_m
Каждый раз проверяем, получился ли 0: gcd(x, m)
================================================
Китайская теорема об остатках
Пусть m = m1 * m2 * m3 * ... * m_k, gcd(m_i, m_j) = 1 при i != j
~
Тогда кольцо Z_m = Z_m1 + Z_m2 + Z_m3 + ... + Z_mk
Прямая сумма двух колец A и B есть множество пар
(a, b)
в котором операции выполняются покомпонентно
(a1, b1) + (a2, b2) = (a1+a2, b1+b2)
Пример: Z_35 = Z_5 + Z_7
8 = (3, 1)
11 = (1, 4) +
-----------
19 = (4, 5)
Алгоритм в Китайской теореме об остатках:
дано k взаимно простых чисел m1, m2, ..., m_k
и k произвольных чисел x1, x2, ..., x_k
Надо найти x: x == x1 (mod m1)
x == x2 (mod m2)
. . .
x == x_k (mod m_k)
Задача: сколько квадратных корней из 1 в кольце Z_105?
x^2 == 1 (mod 105)
105 = 3*5*7
Z_105
x = (x1, x2, x3) x1 <- Z_3, x2 <- Z_5, x3 <- Z_7
1 = (1, 1, 1)
(-1, 1, 1)
(1, -1, 1)
(-1, -1, 1)
. . .
-1 = (-1, -1, -1)
-----------------------------------
Надо найти x: x == x1 (mod m1)
x == x2 (mod m2)
. . .
x == x_k (mod m_k)
Первое равенство
x = x1 + s*m1
Подберем s так, чтобы выполнялось второе равенство
x == x2 (mod m2)
x1 + s*m1 == x2 (mod m2)
s*m1 == (x2 - x1) (mod m2)
Поскольку m1 взаимно просто с m2, то m1 обратимо в кольце Z_m2
Пусть u -- обратный элемент к m1
u*m1 == 1 (mod m2)
Домножим равенство на u:
s*m1*u == (x2 - x1)*u (mod m2) =>
s = (x2 - x1)*u (mod m2)
На шаге i мы вычислили x:
x == x1 (mod m1)
x == x2 (mod m2)
. . .
x == x_i (mod m_i)
Надо подправить x так, чтобы выполнялось еще одно равенство
x' == x_i+1 (mod m_i+1)
Возьмем x' = x + s*(m1*m2*...*m_i)
x + s*(m1*m2*...*m_i) == x_i+1 (mod m_i+1)
Число m1*m2*...*m_i взаимно просто с m_i+1, найдем обратный элемент (mod m_i+1)
u*m1*m2*...*m_i == 1 (mod m_i+1)
s = (x_i+1 - x)*u (mod m_i+1)
=============================================
26 Mar 2022
Классы в Питоне
Питон -- объектно-ориентированный язык. Особенность Питона --
динамическия типизация: типы объектов и выражений невозможно
понять по тексту программы, тип определяется только во время
выполнения; переменные не имеют типов и в разные моменты времени
могут содержать ссылки на объекты разных типов.
Язык вообще не содержит описаний, а только выполняемые операторы.
Классы в Питоне совсем не похожи на классы в других языках.
В C++ класс -- это как бы чертех некоторого объекта (чертеж автомобиля).
Все объекты класса в C++ имеют одно и то же строение (т.е. переменные-члены
класса, объект класса -- это набор переменных-членов + описание
методов, которые работают с объектами класса).
В Питоне всё сосем не так!
Нет никакаго "чертежа" класса, есть только описания методов.
В Питоне разные объекты одного и того же класса могут иметь
разные наборы переменных-членов в зависимости от истории
их создания и модификации. Как правило, переменные-члены создаются
в конструкторе класса, но это не обязательно.
Пример: в C++ класс Вектор на плоскости R2
class R2Vector {
public:
double x;
double y;
public:
R2Vector(double xx = 0., double yy = 0.):
x(xx),
y(yy)
{}
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;
}
double operator*(const R2Vector& v) const { // s = u*v;
return x*v.x + y*v.y;
}
R2Vector operator*(double c) const { // w = v*2.5;
return R2Vectotr(x*c, y*c); // u = 3.33*v; Так написать нельзя!
}
. . .
};
Рассмотрим способы создания класса и работы с классами на примере
класса R2Vector
====================================================
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 __getitem__(self, idx): # x = v[0]
if idx == 0:
return self.x
elif idx == 1:
return self.y
else:
raise IndexError("Incorrect index of vector")
def __setitem__(self, idx, value):
if idx == 0:
self.x = float(value)
elif idx == 1:
self.y = float(value)
else:
raise IndexError("Incorrect index of vector")
def __str__(self):
return "(" + str(self.x) + ", " + str(self.y) + ")"
def __repr__(self):
return "R2Vector" + str(self)
def __add__(self, v):
'''Return self + v'''
return R2Vector(self.x + v[0], self.y + v[1])
def __iadd__(self, v): # in-place addition
'''self += v'''
self.x += v[0]
self.y += v[1]
return self
def dot(self, v): # s = v.dot(w)
'''Dot product self*v'''
return self.x*v[0] + self.y*v[1]
def __mul__(self, v):
'''Return self*v
(either a dot product of 2 vector or a product of vector by a number)'''
if type(v) == R2Vector or type(v) == list or type(v) == tuple:
return self.dot(v)
else:
return R2Vector(self.x*v, self.y*v) # w = v*2.5
# нельзя w = 2.5*v
def __imul__(self, c):
self.x *= c
self.y *= c
return self
def __rmul__(self, v):
'''Return v*self
(either a dot product of 2 vector or a product of vector by a number)'''
if type(v) == R2Vector or type(v) == list or type(v) == tuple:
return v[0]*self.x + v[1]*self.y
else:
return R2Vector(v*self.x, v*self.y) # w = 3.7*u
def __sub__(self, v):
'''Return self - v'''
return R2Vector(self.x - v[0], self.y - v[1])
def __isub__(self, v): # in-place subtraction
'''self -= v'''
self.x -= v[0]
self.y -= v[1]
return self
def __matmul__(self, v):
'''Return a @ v (dot product)'''
return self.dot(v)
def norm2(self):
return self.x*self.x + self.y*self.y
def norm(self):
'''Return the Eucledian norm of vector'''
return self.norm2()**0.5
def __len__(self):
'''Return len(v) = 2'''
return 2
def length(self):
return norm(self)
def normal(self):
'''Return a normal to the vector'''
return R2Vector(-self.y, self.x)
def normalized(self):
l = self.norm()
if l > 0.:
return self*(1./l)
return R2Vector(0., 0.)
def normalize(self):
l = self.norm()
if l > 0.:
self *= 1./l
# return self
====================================================
Решим задачу: найти расстояние от точки до прямой.
Точка: (3, 3)
Прямая через точки (1, 0) и (2, 2)
>>> from R2Vector import *
>>> t = R2Vector(3, 3)
>>> a = R2Vector(1, 0)
>>> b = R2Vector(2, 2)
>>> v = b - a
>>> n = v.normal()
>>> n.normalize()
>>> n
R2Vector(-0.8944271909999159, 0.4472135954999579)
>>> d = abs((t - a)*n)
>>> d
0.44721359549995787
=======================================
2 Apr 2022
Работа с матрицами. Модуль numpy
Матрицы можно представлять как двумерные массивы:
список строк, каждая строка -- список чисел
a = [[1., 2., 3.], [4., 5., 6.], [7., 8., 0.]]
С такими матрицами можно было бы реализовать все матричные
алгоритмы (умножение матриц, приведение к стпенчатому виду,
вычисление обратной матрицы, решение системы линейных уранений, ...).
Однако реализация сложных алгоритмов в Питоне была бы крайне
неэффективной! Решение -- все алгоритмы, все объемные и сложные
вычисления реализуются внутри модулей, которые написаны уже не на Питоне,
а на эффективных языках типа Fortran, C, C++. На Питоне пишется
только самый верхний уровень: создать матрицу, перемножить две матрицы,
вычислить обратную матрицу, решить линейную систему.
Команды отдаются на Питоне, а выполняются уже внутри модуля.
Самый популярный модуль в Питоне -- это модуль numpy
(number computing in python). Я не видел ни одной практической
программы, которая бы не использовала модуль numpy.
Главное, что умеет модуль numpy -- это работать с массивами,
двумерными массивами (матрицами), многомерными массивами и т.д.
import numpy as np
a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 0]], dtype="float64")
Внимание! Важно:
1) все элемента numpy-массива имеют один и тот же тип;
2) типы, если только это не объекты, ограничены набором типов
языка C. Например, целые числа всегда ограничены:
int8, int16, int32, int64;
3) элементы матрицы и вообще многомерного массива хранятся в
линейном массиве:
a00 a01 a02
a10 a11 a12
a20 a21 a22
хранятся по строкам (по умолчанию):
a00 a01 a02 a10 a11 a12 a20 a21 a22
Для numpy-массивов реализована векторная арифметика.
Арифметические операции выполняются покомпонентно.
Скалярное произведение (dot product) выполняется
либо с помощью метода или функции dot
>>> c = np.array([1., 2., 33.])
>>> d = np.array([10., 20., 30.])
>>> c*d
array([ 10., 40., 990.])
>>> c.dot(d)
1040.0
>>> np.dot(c, d)
1040.0
либо с помощью знака @, который в современном Питоне
зарезервирован для скалярного произведения векторов
или произведения матриц:
>>> c.dot(d)
1040.0
>>> np.dot(c, d)
1040.0
>>> c @ d
1040.0
Произведение матриц:
>>> a = np.array([[1, 2], [3, 4]])
>>> a
array([[1, 2],
[3, 4]])
>>> b = np.array([[0, -1], [1, 0]])
>>> b
array([[ 0, -1],
[ 1, 0]])
>>> a.dot(b)
array([[ 2, -1],
[ 4, -3]])
>>> a @ b
array([[ 2, -1],
[ 4, -3]])
Получение миноров матрицы:
>>> a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 0]], dtype="float64")
>>> a
array([[1., 2., 3.],
[4., 5., 6.],
[7., 8., 0.]])
Первая строка матрицы:
>>> a[0]
array([1., 2., 3.])
Минор из первой и последней строки:
>>> a[[0, 2]]
array([[1., 2., 3.],
[7., 8., 0.]])
Второй столбец матрицы:
>>> a[:, [1]]
array([[2.],
[5.],
[8.]])
Элементарные преобразования:
к первой строке прибавить вторую, умноженную на число 10:
>>> a[0] += a[1]*10
>>> a
array([[41., 52., 63.],
[ 4., 5., 6.],
[ 7., 8., 0.]])
Поменять местами две строки матрицы:
>>> a[[0, 2]] = a[[2, 0]]
>>> a
array([[ 7., 8., 0.],
[ 4., 5., 6.],
[41., 52., 63.]])
Умножить строку на число:
>>> a[1] *= -2
>>> a
array([[ 7., 8., 0.],
[ -8., -10., -12.],
[ 41., 52., 63.]])
Это уже достаточно, чтобы реализовать метод Гаусса
================================
9 Apr 2022
Оконное/графическое программирование на Питоне
с использованием модуля tkinter
tkinter -- очень простой оконный интерфейс, буквально
предназначенный для детского сада.
Достоинства -- легко учится, позволяет быстро создавать оконные программы
Недостатки -- не сглаживает линии (antialiasing) (под X-Window;
под Mac OS/X Quartz всё сглаживается);
наверно, очень сложные оконные программы
трудно написать под tkinter.
---------------------------------
16 Apr 2022
Оконное/графическое программирование на Питоне
с использованием модуля tkinter
Рисование графика функции (программа-образец)
. . .
==============================================
23 Apr 2022
Задача линейной регрессии и ее решение с помощью метода
наименьших квадратов. Псевдообратная матрица Мура--Пенроуза
Пусть мы хотим построить линейную зависимость между
независимой переменной x и зависящей от нее переменной y = y(x).
Независимая переменная x <- R^n является точкой n-мерного пространства,
которая описывает некоторый объект. Этот объект характеризуется
набором из n чисел, которые являются признаками объекта (features).
Независимая переменная y <- R -- это просто число.
Зависимость описывается линейной функцией
y = f(x) = <w, x> + b, w <- R^n, b <- R
Пример: хотим построить формулу для стоимости квартиры.
Квартира -- это объект x, который можно описывать следующими произнаками:
1) площадь;
2) количество комнат;
3) этаж;
4) расстояние до ближайшей станции метро;
5) возраст дома;
6) насколько экологичный район;
7) расстояние до центра города;
. . .
x = (x0, x1, ..., x7)
<w, x> = w0*x0 + w1*x1 + ... + w6*x6 + b
w_i -- вес i-го признака.
У нас есть тренировочная выборка
x1, x2, x3, ..., x_m
y1, y2, y3, ..., y_m
y_i -- стоимость i-й квартиры в нашей выборке.
По этой выборке мы хотим построить оптимальную модель,
т.е. найти ее параметры w, b
Почему эта задача называется "регрессией"?
Был ученый, который изучал зависимость роста взрослого сына от роста отца.
x_i -- рост отца, y_i -- рост сына.
Получилась прямая с углом наклона < 45 градусов.
У высокого отца в среднем сын ниже него; и низкого отца, наоборот, выше.
Получается, что с поколениями рост должен усредняться, отсюда и произошло
слово регрессия.
На самом деле, можно упростить зависимость, отказавшись от свободного члена
b. Добавим ко всем объектам дополнительный признак, тождественно равный 1.
x --> x' = (x, 1) <- R^n+1
w --> w' = (w, b)
Зависимость
f(x) = <w, x> + b
переписывается в виде
f(x') = <w', x'>
Слово "линейная" не означает, что мы не можем, допустим,
находить полиномиальную зависимость. Пусть, например,
мы хотим найти зависимость в виде полинома степени d
y = w0*x^d + w1*x^(d-1) + ... + w_d-1*x + w_d
Наш объект описывается степениями переменной x
z = (x^d, x^d-1, ..., x, 1)
y = f(z) = <w, z>
В алгоритмах машинного обучения мы всегда стараемся подобрать
параметры модели так, чтобы на тренировочной выборке она
выдавала бы минимальную ошибку. Как измерять ошибку?
В простейшем случае мы минимизируем среднюю квадратичную ошибку
(Mean Squared Error):
MSE = Sum_i (f(x_i) - y_i)^2
Вариант: можно рассматривать среднюю абсолютную величину ошибки
MAE = Sum_i |f(x_i) - y_i|
Итак, наша задача -- найти вектор w, который минимизирует
средне-квадратичное отклонение.
MSE = Sum_i (<w, x_i> - y_i)^2 --> min_w
Запишем это в матричном виде: координаты x_i записываются в i-й строке
матрицы A, x_i <- R^n, i = 1, 2, ..., m, причем m >= n.
x_11 x_12 x_13 ... x_1n
A = x_21 x_22 x_23 ... x_2n
x_31 x_32 x_33 ... x_3n
. . .
x_m1 x_m2 x_m3 ... x_mn
y_1
Y = y_2
y_3
...
y_m
w_1
W = w_2
w_3
...
w_n
Модель выдает значения A*W. Мы как бы пытаемся решить систему
A W = Y, A -- mxn-матрица, W -- столбец высоты n,
Y -- столбец высоты m >= n
В этой системе m уравнений и n неизвестных w_j, j = 1, 2, ..., n
В общем случае она не имеет решений.
Но тогда попробуем найти "псевдорешение" такое, что расстояние
между векторами AW и Y минимально:
|| A W - Y ||^2 --> min
Это и будет решение задачи регрессии в случае функции ошибки = сумме
квадратов отклонений (или корню квадратному из суммы отклонений).
A -- это линейный оператор
A: R^n --> R^m, m >= n
Надо найти вектор V такой, что расстояние между A V и Y минимально.
То есть вектор A V должен быть проекцией вектора Y на подпространство
Im(A) = A R^n =>
Вектор Y - AV перпендикулярен подпространству Im(A) =>
<(Y - AV), Im(A)> = 0 (скалярное произведение)
(Y - AV)* (A W) = 0 для всякого вектора W
((Y - AV)* A) W = 0 для любого W =>
(Y - AV)* A = 0
(звездочкой обозначаем транспонирование). (S T)* = T* S*
(Y* - V* A*) A = 0
Y* A - V* A* A = =
V* A* A = Y* A
Транспонируем обе части равенства:
A* A V = A* Y
Это называется нормальное уравнение.
Если столбца матрицы A линейно независимы, то матрица A* A обратима
(это так называемая матрица Грамма = матрица попарных скалярных
произведений набора векторов, являющихся столбцами матрицы A).
Тогда мы получаем решение
V = (A* A)^(-1) A* Y
Матрица
A+ = inv(A* A) A*
называется псевдообратной матрицей Мура--Пенроуза. Она обобщает
понятие обратной матрицы для не квадратной матрицы.
A W = Y
квадратная невырожденная матрица A => решение
W = A^(-1) Y
Прямоугольная матрица с числом строк > число столбцов,
столбцы линейно независимы => псевдорешение
W = A+ Y
где A+ -- псевдообратная матрица, A+ = (A* A)^(-1) A*
В Питоне вычисление всех этих матриц реализовано
в модуле numpy.linalg
a_inv = np.linalg.inv(a)
a_pinv = np.linalg.pinv(a)
Матрица Вандермонда:
x1, x2, ..., x_m <- R
Ищем многочлен степени d:
x1^d x1^d-1 ... x1 1
A = x2^d x2^d-1 ... x2 1
x3^d x3^d-1 ... x3 1
. . .
xm^d xm^d-1 ... xm 1
=========================================================
30 Apr 2022
Машинное обучение, элементы искусственного интеллекта
Была рассмотрена задача регрессии
x <- R^n -- независимая переменная x = (x1, x2, ..., x_n)
x_i -- "признаки" объекта x,
features
y = y(x) -- зависимая переменная, y <- R
Мы ищем линейную зависимость
y = <w, x> + b
Имеется тренировочная выборка из m елементов
x1, x2, x3, ..., x_m <- R^n
y1, y2, y3, ..., y_m
Мы подбираем линейную модель, зависящую от параметров w <- R^n, b <- R
так, чтобы она максимально соответствовала нашей выборке:
y = <w, x_i> + b -- значение нашей модели на объекте x_i
y_i -- реальное значение для объекта из тренировочной
выборке.
Ошибка -- это значение y_i - (<w, x_i> + b)
В методе наименьших квадратов мы минимизируем среднюю квадратичную ошибку
на тренировочной выборке:
MSE = 1/m Sum_i (y_i - (<w, x_i> + b))^2 --> min_w,b
MSE -- Mean Squared Error
===============================================================
Рассмотрим вторую типичную задачу машинного обучения:
задачу классификации
Классический метод построения линейного классификатора --
метод опорных векторов (Support Vector Machine).
Пусть мы имеем объекты двух классов, каждый объект описывается
набором из n признаков:
x = (x1, x2, ..., x_n) <- R^n
Мы хотим построить линейный классификатор f(x), который разделяет
эти два класса объектов:
f(x) = <w, x> - b, w <- R^n, b <- R -- свободный член
intercept
Мы считаем, что, если
f(x) > 0, то объект x принадлежит к первому классу,
f(x) < 0, то объект x принадлежит второму классу,
f(x) = 0, то непонятно.
Уравнение
<w, x> - b = 0
определяет гиперплоскость в пространстве R^n коразмерности 1.
Эта гиперплоскость должна разделять наши 2 класса точек:
точки в верхнем полупространстве принадлежат 1-му классу,
в нижнем -- второму классу.
Как построить такую гиперплоскость?
Рассмотрим сначал разделимый случай, когда такая гиперплоскость
существует. Предполагаем, что выполнена СИММЕТРИЧНАЯ разметка
данных:
x1, x2, x3, ..., x_m
y1, y2, y3, ..., y_m = +-1
y_i = { 1, если x_i принадлежит 1-му классу
{ -1, если x_i принадлежит 2-му классу
Тогда, если классификатор выдает правильный ответ на объекте x_i, то
величина положительна:
y_i (<w, x_i> - b) > 0
Утверждение: расстояние со знаком от точки x до гиперплоскости,
заданной уравнением
<w, x> - b = 0,
выражается формулой
h = (<w, x> - b)/|w|
Почему это так? w -- это вектор нормали к гиперплоскости, но длина его
не обязательно равна 1. Тогда v = w/|w| -- единичная нормаль и
<w/|w|, x> - b/|w| = h
(<w, x> - b)/|w| = h
---------
Для опорных векторов x_i имеем точное равенство
y_i(<w, x_i> - b)/|w| = h
Для остальных векторов имеем неравенство
y_j(<w, x_j> - b)/|w| > h
То есть для любых векторов должно быть выполнено условие
y_i(<w, x_i> - b)/|w| >= h для всех i=1, 2, ..., m
Мы хотим также максимизировать ширину полосы, разделяющей
точки двух классов: ширина полосы равна 2*h
Преобразуем неравенство
y_i(<w, x_i> - b)/|w| >= h
Разделим обе его части на h
y_i(<w/(h*|w|), x_i> - b/(|w|*h)) >= 1
Обозначим
w' = w/(h*|w|), b' = b/(h*|w|)
Получим неравенство
y_i(<w', x_i> - b') >= 1
Выразим h через |w'|:
w' = w/(h*|w|)
|w'| = |w|/(h*|w|)
|w'| = 1/h
h = 1/|w'|
Поучили, что ширина полосы h = 1/|w'|
Нам надо максимизировать h
Получаем задачу (в обозначениях опустим штрихи):
y_i(<w, x_i> - b) >= 1, i=1, 2, ..., m
|w|^2 --> min_w,b
Это задача метода опорных векторов в разделимом случае:
минимизировать функцию |w|^2 при заданном наборе ограничений
в виде неравенств
y_i(<w, x_i> - b) >= 1, i=1, 2, ..., m
Это так называемая задача квадратичного программирования.
Она решается не всегда, потому что не всегда можно разделить точки
двух классов гиперплоскостью.
Как поступить???
Давайте ослабим неравенства
y_i(<w, x_i> - b) >= 1, i=1, 2, ..., m
введя дополнительные переменные (расслабляющие переменные,
slack variables) xi_i, i = 1,..., m
y_i(<w, x_i> - b) >= 1 - xi_i
xi_i >= 0
Число xi_i -- это как бы величина штрафа за то, что мы не можем точно
выполнить i-е неравенство.
Мы минимизируем штрафную функцию
|w|^2 + C*(1/m)*Sum xi_i --> min
при условиях
y_i(<w, x_i> - b) >= 1 - xi_i,
xi_i >= 0, i=1,2,...,m
======================================
7 May 2022
Нелинейный вариант метода опорных векторов
Мы рассматривали задачу регрессии
Есть независимая переменная x <- R^n x = (x1, x2, ..., x_n)
x_i <- R -- признаки объекта
(features)
В задаче регрессии есть зависимая переменная y = y(x)
причем зависимость линейная:
y = <w, x> + b w <- R^n -- вектор весов
Есть тренировочная выборка
x1, x2, ..., x_m <- R^n
на этих объектах функция принимает значения
y1, y2, ..., y_m <- R
Параметры модел w, b подбираются так, чтобы значения, которая она выдает
на векторах x_i, были бы максимально близки к y_i.
Было рассмотрено, как эта задача решается методом наименьших квадратов
с помощью псевдообратной матрицы.
Задача классификации:
есть два класса объектов x <- R^n. Мы хотим построить линейный классификатор
f(x) = <w, x> - b
так, что он относит объект x к паервому классу, когда f(x) > 0,
и ко второму классу, кгда f(x) < 0.
Когда f(x) = 0, то непонятно, к какому классу отнести x.
Пусть мы имеем тренирововчную выборку
x1, x2, ..., x_m <- R^n
причем известно, к какому классу принадлежит каждый из этих объектов:
y1, y2, ..., y_m = +-1
{ 1, если x_i принадлежит первому классу
y_i = {
{ -1, если x_i принадлежит второму классу
Уравнение
<w, x> - b = 0
задает гиперплоскость в n-мерном пространстве, разделяющую эти
два класса точек. Эта гиперплоскость в случае, когда это сделать можно,
находится как решение оптимизационной задачи
y_i(<w, x_i> - b) >= 1 для всех i=1,2,...,m
|w|^2 --> min
Если разделить классы гиперплоскостью нелья, то гиперплоскость находится
как решение оптимизационной задачи:
y_i(<w, x_i> - b) >= 1 - ksi_i для всех i=1,2,...,m
ksi_i >= 0
|w|^2 + (C/m) * Sum ksi_i --> min
Из этих неравеств можно исключить переменные ksi_i:
y_i(<w, x_i> - b) >= 1 - ksi_i для всех i=1,2,...,m
ksi_i >= 0
=>
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)
Подставив начение ksi_i в целевую функцию, получим задачу безусловной
минимизации:
|w|^2 + (C/m) * Sum max(1 - y_i(<w, x_i> - b), 0) --> min
Используем функцию Hinge Loss
H(c) = max(1-c, 0)
Получаем задачу:
|w|^2 + (C/m) * Sum H(y_i(<w, x_i> - b)) --> min_w,b
В прошлый раз мы ее реализовали на Питоне.
Как быть, когда точки двух классов невозможно разделить
гиперплоскостью?
Идея: давайте дополним признаки объектов некоторым набором
базисных функций от исходных признаков:
вместо x = (x1, x2, ..., x_n) возьмем
x' = (phi_1(x), phi_2(x), ..., phi_k(x)), k >= n
Отображение phi переводит точку x в пространство большей размерности:
phi: x --> x'
phi: R^n --> R^k, k >= n
Можно надеяться, что в пространстве большей размерности точки
уже можно разделить гиперплоскостью.
F(x') = <W, x'> - B
Классификатор для исходного пространства выражается как
f(x) = F(phi(x)) = <W, phi(x)> - B
Пример.
В качестве базисных функций можно взять все мономы (произведения)
степени <= d от исходных признаков. Например, пусть объект
p = (x, y)
тогда
p' = (1, x, y, x^2, x*y, y^2)
phi: R^2 --> R^6
Если степень d = 3, то
p' = (1, x, y, x^2, x*y, y^2, x^2, x^2*y, x*y^2, y^3)
Как получить все мономы?
Воспольуемся модулем itertools и в нем генераторной функцией
combinations_with_replacement(x, d)
которая генерирует все упорядоченные наборы длины d
от элементов списка x (воможно, с повторениями)
>>> import numpy as np
>>> import itertools
>>> x = [2, 3, 5]
>>> for p in itertools.combinations_with_replacement(x, 2):
... print(p)
...
(2, 2)
(2, 3)
(2, 5)
(3, 3)
(3, 5)
(5, 5)
>>> for p in itertools.combinations_with_replacement(x, 2):
... print(p)
... print(np.prod(p))
...
(2, 2)
4
(2, 3)
6
(2, 5)
10
(3, 3)
9
(3, 5)
15
(5, 5)
25
Все мономы степени <= d генерируются следующей функцией:
def monomials(x, degree=2):
res = [1.]
for d in range(1, degree + 1):
for p in itertools.combinations_with_replacement(x, d):
res.append(np.prod(p))
return res