Язык Crystal и математика

Если вам вдруг захотелось посчитать собственные значения матрицы, решить задачу линейного программирования или оптимизировать нелинейную функцию, то вы может взять Python со SciPy, можете взять R или Matlab/Octave, для любителей экзотики есть Julia, а те кому важен каждый тик скорее возьмут C++ или Rust.

Но если вдруг вы захотите взять Crystal (оставим за кадром причины этого решения), то я постараюсь описать что он может предложить для решения математических задач.

1. Стандартная библиотека

2. Использование библиотек с С API

3. Мои библиотеки

   1. linalg

   2. crystal-gsl

   3. nlopt

   4. linprog

4. Не мои библиотеки

   1. num.cr

   2. fftw.cr

Стандартная библиотека

Crystal унаследовал от Ruby приличное количество "батареек" в стандартной библиотеке. Вкратце пробежимся по тому что можно отнести к математике.

Модуль Math включает в себя 41 математическую функцию - почти все вызывают libm или интринсики llvm.

Модуль Complex позволяет не только складывать комплексные числа, но и посчитать от них корень, экспоненту и логарифм. Внутри состоят из двух Float64.

BigInt / BigFloat / BigDecimal / BigRational позволяют работать с числами с неограниченной размерностью (внутри используют известную библиотеку GMP).

В модуле Indexable есть такие чудеса как each_permutation и даже #each_repeated_combination - не то чтобы их сложно было написать самому, но приятно что их уже написали за тебя. А, ну еще #tally не могу не упомянуть - мелочь, но удобная.

Использование библиотек с С API

В Crystal легко создавать байндинги к библиотекам с С API. Это конечно не C/C++ где всё уже готово и надо только слинковаться и не Go, где cgo позволяет импортировать библиотеку парой строк, но и не Ruby, где при подключении сишной библиотеки вся надежда на то что кто-нибудь уже написал соответствующий gem, не то придется колдовать с Fiddle и обертками к каждой функции.

Crystal в этом смысле больше всего похож (из знакомых мне языков) на паскаль - тоже генерируем заголовочный файл по h файлу, проходимся ручками чтоб исправить косяки и вуаля - можно звать LibEngine.ellipse(center.x, center.y, radius.x, radius.y, filled, color1, color2, angle) из своего кода.

Конкретно генерация выполняется библиотекой crystal_lib. Внутри она использует clang для синтаксического анализа. Преимущество (по сравнению с парсером на регулярных выражениях) - хитрые синтаксические выверты для нее не проблема, недостаток - для получения результата заголовок должен успешно скомпилироваться.

Что это означает для математического применения? Что если есть интересующая математическая библиотека (обязательно с C, не С++ API), не так уж сложно использовать ее из Кристалла.

Мои библиотеки

Linalg

Моя библиотека для работы с матрицами, вдохновленная matlabом и scipy. Точнее так - сначала я реализовал в ней всё чем пользовался в свое время в матлабе, а дальше ориентировался на список функций scipy.linalg.

Ссылка - https://github.com/konovod/linalg/

неполная автоматически сгенеренная документация https://konovod.github.io/linalg/

Сразу важное замечание. Говоря о матрицах в программировании подразумевают два вида библиотек. Один - для матриц 4*4 используемых в компьютерной графике, второй - для матриц размером 1000*1000. Формально конечно операции с ними одни и те же, но акценты при разработке разные. Я интересуюсь вторым типом.

Как реализована библиотека - простые вещи типа triu реализованы на кристалле, умножение матриц выполняется с помощью BLAS, а сложная линейная алгебра - с помощью LAPACK. На винде и линуксе для ее работы потребуется поставить OpenBLAS, на маке используется встроенная в систему реализация (фреймворк Accelerate).

Для тех кто не знает что такое.

BLAS - стандартный интерфейс библиотек для умножения матриц (плотных, т.е. dense). Есть разные реализации, из опенсорсных самая "продвинутая" по моим данным OpenBLAS, из проприетарных есть Intel MKL, NVidia cuBLAS, Apple Accelerate.

LAPACK - это внушительная коллекция процедур линейной алгебры на фортране. Примерно 5 сотен процедур (в последней версии больше 7 сотен), каждая из которых доступна в 2 или 4 вариантах (для двойной и одинарной точности, для действительных и комплексных чисел). Все называются непонятными для непосвященных аббревиатурами (например dgees, csytrf_aa или zggbal). Все подробно документированы и тщательно протестированы. Дополняется и исправляется по сей день (https://github.com/Reference-LAPACK/lapack). Большинство из этих 5 сотен служебные (впрочем, строго провести границу между служебными и неслужебными не так легко), некоторые устарели, но и оставшихся хватает. В моей библиотеке я использую 57 функций оттуда, еще несколько десятков ждут своей очереди.

Возвращаясь к моей библиотеке - что можно делать с ее помощью.

• создание матриц, доступ к элементам по индексу и диапазону, map и tril

• форматированный ввод\вывод

• подматрицы

• арифметика

• линейная алгебра - det, inv, solve, norm, rank

• недо и переопределенные задачи

• собственные значения (eigenvalues), SVD, псевдоинверсия

• декомпозиция матриц.(LU, QR, RQ, LQ, QL, cholesky, schur\QZ)

• флаги матриц

  Флаги - в каждой матрице хранятся специальные свойства которыми она обладает (треугольная, симметричная, эрмитова, положительно определенная, ортогональная), эти свойства учитываются при операциях (например треугольную можно умножать быстрее, а определитель ортогональной всегда 1).

• ленточные матрицы (banded matrix)

  Это частный случай разреженных - матрицы у которых хранится только главная диагональ и н диагоналей выше и ниже нее. Не знаю насколько они распостранены, но в LAPACK им посвящено приличное количество функций, так что я добавил их поддержку.

• специальные матрицы

• матричная экспонента.

• разреженные матрицы (sparse matrix)

  С ними все неидеально. LAPACK все-таки ориентирован на плотные матрицы, какую библиотеку выбрать для разреженных я пока не решил (все оборачивать жизни не хватит, на бесполезные время тратить не хочется, пока на примете SuiteSparse но там нет svd и некоторых других функций). Пока на чистом кристалле сделал хранение и преобразование матриц COO, CSR.

crystal-gsl

Ссылка - https://github.com/konovod/crystal-gsl

автоматически сгенеренная документация https://konovod.github.io/crystal-gsl/

В отличие от прошлой библиотеки, где кроме вызовов LAPACK пришлось написать кучу своего кода, эта - просто обертка над GSL. Ну, какой-то код там естественно пришлось добавить, но большая часть работы была в том чтобы придумать красивый ООП интерфейс к каждой функции и написать тесты.

Что такое GSL (https://www.gnu.org/software/gsl/doc/html/index.html). Библиотека для математических расчетов от проекта GNU. Распостраняется под лицензией GPL. Несмотря на внушительный список областей, общее впечатление довольно грустное. Даже LAPACK который пишут на фортране мигрировал на гитхаб и динамично развивается, а у GSL один активный разработчик и баги с уже приложенным патчем для исправления висят в трекере годами и десятилетиями. Да-да, опенсорс убил спо. Но хватит о грустном, что-то хорошее там все-таки есть, так что если мы хотим конкурировать с scipy, берем, пишем обертки и добавляем сахар.

Изначально обертку сделал ruivieira (https://github.com/ruivieira/crystal-gsl), но он потерял интерес к проекту (а я наоборот приобрел), так что актуальный форк теперь мой.

Что там есть

• сгенерирован файл для актуальной версии GSL, так что можно вызывать все функции напрямую.

• ООП интерфейс для следующих тем:

  • Общие

    • Векторы и Матрицы

      Для матриц у меня уже есть linalg, так что тут я особого сахара не добавлял, но те функции какие в GSL есть обернуты, складывать и умножать матрицы можно.

      Пример

      m1 = GSL::DenseMatrix.new(3, 3)
      m1[0, 2] = 2.0
      m2 = [  [1.0, 0.6, 0.0],
              [0.0, 1.5, 1.0],
              [0.0, 1.0, 1.0],
            ].to_matrix
      puts (m1*m2.inverse)

    • Разреженные (sparse) и Плотные(dense) матрицы.

      Вот разреженные матрицы - то чего мне в linalg не хватает. В GSL уровень конечно сильно отстает от современных достижений (типа SparseSuite), но это хоть что-то! Собственно, что там поддерживается: типы хранения COO, CSC, CSR, конверсия между ними, перемножение матриц (только CSC на CSC) и единственный итеративный солвер GMRES.

      Пример

      describe ".solve" do
        it "should solve example from gsl docs" do
          n = 100         #  number of grid points
          size = n - 2    # subtract 2 to exclude boundaries
          h = 1 / (n - 1) # grid spacing
          a = GSL::SparseMatrix.new(size, size)
          f = GSL::Vector.new(size)
          # construct the sparse matrix for the finite difference equation
          # construct first row
          a[0, 0] = -2
          a[0, 1] = 1
          # construct rows [1:n-2]
          (1...size - 1).each do |i|
            a[i, i + 1] = 1
            a[i, i] = -2
            a[i, i - 1] = 1
          end
          # construct last row
          a[size - 1, size - 1] = -2
          a[size - 1, size - 2] = 1
          # scale by h^2
          a = a*(1/h/h)
          # construct right hand side vector
          size.times do |i|
            xi = (i + 1) * h
            fi = -Math::PI * Math::PI * Math.sin(Math::PI * xi)
            f[i] = fi
          end
          # convert to compressed column format
          c = a.convert(:csc)
          # now solve the system with the GMRES iterative solver
          u = GSL::SparseMatrix.solve(a, f)
          # compare with exact solution
          size.times do |i|
            xi = (i + 1)*h
            exact = Math.sin(Math::PI*xi)
            u[i].should be_close exact, 1e-3
          end
        end
      end

  • Статистика. Это не совсем моя область интересов, но этим занимался изначальный создатель байндингов.

    • Гистограммы (не рисуем, только формируем массивы)

    • Перестановки (Permutations) Правда как я упоминал в начале это уже есть в стандартной библиотеке Кристалла

    • Корреляция (Correlation) Пирсон и Спирмен, что бы это ни значило

    Я продолжил работу обернув с помощью макросов

    • 33 из 38 Статистических распределений (distributions)

  • Мат. анализ

    • Численное интегрирование (Numerical integration) более-менее простые алгоритмы, зато протестированы. вызываем GSL.integrate(myfunction, a, b, epsabs : 1e-6) и готово. Можно выбрать алгоритм, некоторые подходят для несобственных интегралов.

    • Численное дифференцирование (Numerical differentiation) GSL.diff(myfunction, x, 0.01) и получаем производную в x по четырем точкам (x-0.01, x-0.005, x+0.005, x+0.01). Есть вариант который берет точки только слева или только справа от x: GSL.diff(myfunction, x, 0.01, :forward).

    • ОДУ (Ordinary differential equations) Есть алгоритмы требующие якобиана и не требующие. Пример: https://github.com/konovod/crystal-gsl/blob/5a79d65f0464f7f095ada9828d51eb1c5b54d863/spec/ode_spec.cr

    • Полиномы Можно хранить полиномы в виде списка коэффициентов, брать от них интеграл и производную и конечно искать корни. Еще есть PolyDD - форма представления полинома в виде Разделённой разности (divided-difference representation)

    • Monte Carlo Integration Интегрирование многомерных функций с помощью метода Монте-Карло. Задаем функцию, границы, получаем приближенный результат. Есть алгоритмы MISER и VEGAS, подробнее о них можно почитать в документации GSL: https://www.gnu.org/software/gsl/doc/html/montecarlo.html

    • Series Acceleration Некий Levin u-transform, позволяющий оценить сумму ряда по первым нескольким членам.

  • Оптимизация

    • Минимизация функций одной переменной. После ознакомления с этим разделом мое мнение о GSL сильно упало. Мало того что вершиной технологий в нем является метод золотого сечения, так он еще и реализован с ошибкой. Баг-репорт об этой ошибке висит с 2015 года: http://savannah.gnu.org/bugs/?45053, так что очевидно что никто этим разделом не пользуется. Эх, вот бы какой-нибудь энтузиаст забил на все эти мэйллисты и форкнул GSL на гитхабе - я уверен что разработка бы сразу в гору пошла. Но тем не менее, я исправил эту ошибку используя приложенный патч и в crystal-gsl таки можно искать минимум функции одной переменной.

    • Поиск корней функции одной переменной Если кому-то интересно чем это отличается от минимизации - при поиске корней у нас в начале интервала у функции один знак, а в конце - другой. Так что уж один корень точно есть и мы его найдем.

    • Алгоритм имитации отжига (Simulated annealing) метод оптимизации, заключающийся в том что мы сначала даем x случайные приращения и ищем минимум, а потом постепенно уменьшаем разброс этих приращений. Реализация состоит из короткой функции, поэтому конкретно GSL тут не используется - я просто переписал эту функцию на Crystal.

    • Линейная регрессия В GSL есть довольно впечатляющее количество функций линейной и нелинейной регрессии, но я пока реализовал только обертку для линейной регрессии функции одной переменной.

  • Аппроксимация

    • Полиномы Чебышева создаем объект poly = GSL::Chebyshev.new(myfunction, a, b) и дальше можем вычислять его в любой точке poly.eval((a+b)/2)также оценивать погрешность и находить производную и интеграл

    • B Сплайны Создаем объект передавая ему значения на равномерной сетке и дальше можем получать его значения и производные в произвольной точке. Эту часть несколько улучшили в 2.8, но текущая стабильная версия GSL 2.7, так что улучшения я пока не оборачивал.

    • Интерполяция идеологически отличается от прошлого пункта тем что интерполирующая функция точно проходит через все заданные точки, а B сплайн лучше подходит когда исходные данные зашумлены. Поддерживается 1- и 2-мерная интерполяция.

  • Физические константы Значения слегка устарели, зато есть в двух системах (МКСА и СГС)

  • Специальные функции 273 из 338 специальных функций обернуты чтобы их можно было вызывать не парясь с синтаксическим мусором. По ощущениям, это самая популярная часть библиотеки GSL - я видел байндинги к ней, которые ограничивались этой частью. Тот же octave, как я понимаю, использует из GSL только эту часть, специальные функции.

• Написаны в папке spec тесты ко всем этим областям. Их же можно использовать и как примеры использования.

  таблица со статусом того что уже обернуто\еще нет\не планируется: https://github.com/konovod/crystal-gsl/blob/master/TODO.md

  На этом с обзором библиотеки GSL закончим, и перейдем к

nlopt

Обертка над библиотекой нелинейной оптимизации NLopt (https://nlopt.readthedocs.io/en/latest/). Библиотека позволяет решать задачи оптимизации с нелинейной целевой функцией и ограничениями в виде линейных и нелинейных неравенств.

Ссылка: https://github.com/konovod/nlopt.cr

Документация: да в общем-то в ридми и приведена.

# создаем объект-решатель
s1 = NLopt::Solver.new(NLopt::Algorithm::LdMma, 2)
# устанавливаем требуемую точность по x (есть еще по f)
s1.xtol_rel = 1e-4
# задаем целевую функцию
s1.objective = ->(x : Slice(Float64), grad : Slice(Float64)?) do
  if grad
    grad[0] = 0.0
    grad[1] = 0.5 / Math.sqrt(x[1])
  end
  Math.sqrt(x[1])
end
a = 2
b = 0
# задаем ограничение в форме неравенства
s1.constraints << NLopt.inequality do |x, grad|
  if grad
    grad[0] = 3 * a * (a*x[0] + b) * (a*x[0] + b)
    grad[1] = -1.0
  end
  ((a*x[0] + b) * (a*x[0] + b) * (a*x[0] + b) - x[1])
end
# задаем ограничение в форме равенства
s1.constraints << NLopt.equality do |x, grad|
  if grad
    grad[0] = 1.0
    grad[1] = 1.0
  end
  x[0] + x[1] - 3
end
# а еще можно задать жесткие границы для переменных и начальное приближение
s1.variables[0].set(min: 0.0, guess: 2.5)
# решаем уже
res, x, f = s1.solve
# res будет NLopt::Result::XtolReached если мы достигли точности по x
# x - вектор значений найденного минимума
# f - значение функции в найденном минимуме

В целом эта библиотека, наверное, не подойдет для супербольших задач, но для любительского использования и не очень больших работает хорошо.

Раньше я немного занимался оптимизацией задач побольше (таких где ответа надо было ждать часами и использовать дорогущий коммерческий солвер, т.к. с опенсорсным ждать пришлось бы днями), ну и вынес оттуда, что лучше всего мы умеем решать линейные задачи, а нелинейные часто можно свести к ним. Как минимум, линейная задача может дать начальную точку для поиска нелинейной. Поэтому следующая библиотека -

Linprog

Это моя обертка над библиотеками для решения ЗЛП (задач линейного программирования).

Ссылка: https://github.com/konovod/linprog/

я поискал какая библиотека самая быстрая и какую при этом легко обернуть и выбрал SYMPHONY (https://projects.coin-or.org/SYMPHONY). У них там кстати на сайте целый мир проектов по оптимизации, причем живые уже перенесли разработку на гитхаб, что радует.

Как выяснилось, SYMPHONY может и быстрая, но сделана специально для "серьезных" задач, так что для мелких задач большую часть времени будет занимать запуск библиотеки. Она еще и пишет в консоль ход решения (в последней версии дали возможность это отключать) - это и неудивительно, если учесть что она ориентирована на задачи которые решаются часами.

В качестве простой альтернативы Linprog я решил добавить GLPK от проекта GNU. Эта библиотека как и GSL в полумертвом состоянии, последние обновления много лет назад, и не потому что нечего улучшать, вместо этого в документации мы с грустью читаем

Note that the GLPK branch-and-cut solver is not perfect, so it is unable to solve hard or very large scale MIP instances for a reasonable time

Но хватит о плохом - зато ее легко установить и она быстро запускается. То что надо для маленьких задач то что надо. Плюс собрать SYMPHONY под винду целый квест, а GLPK собирается без проблем.

В целом linprog пока в не столь развитом состоянии как две мои предыдущие библиотеки - даже CI не настроен. Зато я сделал там DSL для задач, потому что задавать ЗЛП в виде матриц это ад.

# Создаем проблему
task = LinProg::SymbolProblem.new
# Добавляем переменные. параметром можно сразу передать ограничения 
# на целочисленность и мин\макс
x = task.create_var(bound: LinProg::Bound.integer)
y = task.create_var(bound: LinProg::Bound.new(0.0, 6.0, true))
# Добавляем неравенства 
task.st(x + 1 >= y)
task.st(3*x + 2*y <= 12)
task.st(2*x + 3*y <= 12)
# Устанавливаем целевую функцию. Можно и minimize
task.maximize(y + 1)
# Решаем
task.solve
# Теперь из переменных можно прочитать значения
x.value.should eq 1
y.value.should eq 2

Не мои библиотеки

Чтобы не сложилось впечатление что пост только о моих библиотеках, напишу что еще знаю интересного.

Про нейросети писать не буду - я в них недостаточно разбираюсь, да и не особо интересуюсь. Вот бессмысленно и беспощадно оборачивать библиотеку которая ищет SVD для разреженных матриц - это да, это интересно, а всей этой генерацией текста и распознаванием картинок - без меня есть кому заняться.

num.cr

ссылка: https://crystal-data.github.io/num.cr/

Интересный проект который вводит многомерный Tensor и позволяет работать с ним на gpu (с помощью OpenCL). Мне не очень нравится "динамичность" его тензоров - то что даже число измерений задается в рантайме, мой "многомерный массив мечты" выглядит иначе, но то что эта библиотека шире чем моя linalg - несомненно. Кратко о ее возможностях:

• многомерные массивы чего угодно (Tensor)

• могут храниться на цпу и на гпу

• поддерживают эффективный map для расчетов без создания промежуточных объектов

• поддерживают эйнштейновскую нотацию для еще более эффективного вычисления

• линейная алгебра (использует LAPACK)

• градиенты и нейросети

fftw.cr

ссылка: https://github.com/firejox/fftw.cr

Обертка для библиотеки Fastest Fourier Transform in the West (https://fftw.org/)

Я, к моему сожалению, про fft почти ничего не знаю (ну кроме того что это очень круто), но т.к. оно в списке того что есть в GSL - сразу скажу что то что в gsl полная хрень, используйте вот эту быструю и надежную библиотеку. В обертке не хватает некоторых продвинутых возможностей fftw, но т.к. я слабо в этом разбираюсь то не берусь судить насколько они важны.

require "fftw.cr"

x = Array.new(512) { Complex.new(Random.next_u, Random.next_u) }

dft_x = FFTW.dft(x)

require "fftw.cr"

plan = FFTW::Plan.new(512)

x = Array.new(512) { Complex.new(Random.next_u, Random.next_u) }

dft_x = plan.dft(x)