Язык Crystal и математика13.12.2023 17:00

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

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

Стандартная библиотека
Использование библиотек с С API
Мои библиотеки
1. linalg
2. crystal-gsl
3. nlopt
4. linprog
Не мои библиотеки
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) Пирсон и Спирмен, что бы это ни значило
  Я продолжил работу обернув с помощью макросов
- Мат. анализ
  - Численное интегрирование (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 для задач, потому что задавать ЗЛП в виде матриц это ад.

Зачем вам AMPL когда вы можете задавать задачи сразу на любимом языке!

# Создаем проблему
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)

.