Более новая версия находится по адресу: http://wordsandbuttons.online/how_much_math_can_you_do_in
Alines_of_python.html.
Сначала небольшой отказ от ответственности. Этот пост не о линейной алгебре в Python. Если вы хотите именно этого, взгляните на numpy.linalg, оно того стоит. Этот пост о том, насколько прекрасен Python по своей сути и как многого вы можете достичь, написав очень мало кода.
Мы начнем медленно. В Python есть понимание списков. Это способ превратить один список в другой без явного указания циклов и индексов. Итак, если вам нужна функция, которая умножает вектор на скаляр, вы можете сделать это так просто:
def scale(A, x): return [ai*x for ai in A]
Вы можете сделать несколько списков одновременно с помощью zip. Он берет пару списков и делает из них один список кортежей. Тогда каждый кортеж содержит ровно один элемент из каждого списка. Таким образом, добавление векторов тоже довольно просто.
def add(A, B): return [ai+bi for (ai, bi) in zip(A, B)]
Python также имеет некоторые функции обработки списков из коробки. Как 2_. Давайте используем его для написания функции для продукта dot.
def dot(A, B): return sum([ai*bi for (ai, bi) in zip(A, B)])
Понимание списков — это здорово, но вы все равно можете время от времени работать с индексами. В Python есть функция range, которая генерирует для вас последовательность индексов. В Python также есть условные выражения, поэтому с помощью этих двух функций вы можете составить матрицу идентичности в виде вложенного списка.
def ident(n): return [[1.0 if i==j else 0.0 for j in range(n)] for i in range(n)]
Помните zip? Мы использовали его для добавления двух векторов, но он может работать с произвольным количеством массивов. Матрица, как вложенный список, представляет собой произвольный список строк матрицы. Если нам удастся его заархивировать, мы получим список столбцов. Это матричная транспозиция. Кроме того, Python позволяет превратить список в список аргументов для функции с помощью специального синтаксиса. Вы должны написать свой список со звездочкой перед ним.
def transpose(A): return [list(aj) for aj in zip(*A)]
Конечно, вы можете использовать свои собственные функции для создания новых. До сих пор все наши функции были ссылочно прозрачными, что означает, что они всегда обеспечивают один и тот же вывод для одного и того же ввода. Как математические функции или точно настроенные механические машины. Это отличная черта, позволяющая компоновать код без особого усложнения. Мы можем повторно использовать функцию dot ранее в контексте матриц путем умножения векторов.
def mul(A, X): return [dot(ai, X) for ai in A]
Тем не менее, это может стать грязным. В какой-то момент код состоит из определяемых пользователем слов больше, чем из базового синтаксиса. Добавляя новые функции, вы создаете свой собственный язык, поэтому именование так важно.
Мы уже определили 3 слова для векторных операций, давайте объединим их, чтобы создать функцию «проект». Он проецирует точку из пространства на плоскость, гиперплоскость или линию, определяемую ее вектором нормали и расстоянием до центра координат.
def project(A, Pn, Pd): return add(A, scale(Pn, (Pd — dot(Pn, A)) / dot(Pn, Pn)))
Python также предоставляет синтаксис для умножения списков путем повторяющейся конкатенации. Это когда вы пишете [1, 2, 3] * 3 и получаете [1, 2, 3, 1, 2, 3, 1, 2, 3]. Мы можем использовать это для определения функции distance, поскольку скалярное произведение вектора на тот же самый вектор представляет собой просто квадрат его длины.
def distance(A, B): return pow(dot( *(add(A, scale(B, -1.0)), )*2 ), 0.5)
Теперь все может стать немного сложнее. Я хочу показать вам, как решить систему линейных уравнений с. Это не мелочь, у меня ведь целая лекция на эту тему, а люди поумнее меня написали целые книги. Что еще хуже, мы хотим сделать это в одной строке кода.
Это было бы довольно просто с правилом Крамерса, но вы столкнулись бы с серьезными вычислительными проблемами, используя его без конденсации Хио, и это может потребовать поворота, и это не имеет значения для одной линии. У меня есть другая идея.
Помните, вы можете определить плоскость с помощью вектора нормали и скаляра? Ну, это в значительной степени утверждает линейное уравнение. Вы можете представить систему линейных уравнений в виде пересекающихся плоскостей. Когда две из них параллельны, нет ни пересечения, ни решения. Когда два из них идентичны, у вас есть целая линия в качестве решения, поэтому нет особой точки. Но когда они пересекаются, есть одна единственная точка, и вы хотите ее.
И эта точка, очевидно, лежит сразу на каждой плоскости. А теперь подумайте: точка в пространстве не может быть ближе к точке на плоскости, чем ее собственная проекция на эту плоскость. Потому что точка в пространстве, ее проекция и точка на плоскости составляют прямоугольный треугольник, а гипотенуза не может быть короче катета. Поэтому, если мы начнем с любой точки пространства и многократно проецируем ее на каждую плоскость из системы, мы в конце концов приблизимся к решению. Это если он есть конечно.
Так устроены все итерационные алгоритмы. Вы следите за тем, чтобы каждая операция приближала вас к цели, которая называется конвергенцией, и выполняете эти операции до тех пор, пока не закончите.
Но мы хотим, чтобы это было однострочником, поэтому циклы здесь неуместны. Мы будем использовать рекурсию. На каждом шаге мы будем проверять, хорошо ли мы справляемся с текущим решением, а если нет — поворачиваем нашу систему на одно уравнение и идем дальше.
Что касается поворота матрицы, Python предоставляет очень удобный синтаксис для диапазонов списков, который мы с удовольствием воспользуемся. Вы можете не только получить элемент из списка по его индексу в Python, но также указать диапазон элементов, которые вы хотите получить, с помощью двоеточия. Таким образом, все элементы A, кроме первого, будут просто A[1:] , а список одного первого элемента A будет [A[0]]. Теперь вы просто объединяете их, как A[1:] + [A[0]], и вы получаете ротацию.
Это все, что вам нужно для решения системы:
def solve(A, B, Xi): return Xi if distance(mul(A, Xi), B) < 0.0001 else solve(A[1:] + [A[0]], B[1:] + [B[0]], project(Xi, A[0], B[0]))
Давайте завершим наше путешествие вычислением обратной матрицы. Я покажу вам хороший трюк для этого.
Обычно вам нужно найти присоединенную матрицу кофакторов и разделить ее на определитель, но это слишком много печатать. Мы также можем найти обратную матрицу с однострочной.
Учти это. Когда мы умножаем матрицу на вектор, у нас есть список скалярных произведений. Теперь, если вектор является ортом, например [1, 0, 0, …], умножив на него матрицу, мы получим список его элементов первого столбца. Что является первым столбцом обратной матрицы. Мы можем получить любой столбец, выбрав соответствующий орт.
Мы можем получить вектор; мы бы получили, умножая обратную матрицу на что-то, решая систему линейных уравнений, состоящую из исходной матрицы, которая сама является обратной к обратной.
Поэтому мы просто решим эту систему для каждого орта: [1, 0, 0, …], [0, 1, 0, …], [0, 0, 1, …] и так далее, а затем транспонируем список результирующих решений, чтобы получить обратную матрицу.
def invert(A): return transpose([solve(A, ort, [0.0]*len(A)) for ort in ident(len(A))])
Вот и все. Всего 10 маленьких функций, но столько вычислений!
Конечно, повседневный код Python не так уж и лаконичен. Вы можете в значительной степени писать на Python в стиле Fortran или Java, и под «вы», конечно же, я подразумеваю всех остальных. Но что касается основного языка, Python, безусловно, может предложить много интересного.
