Поговорим об отложенных(ленивых) вычислениях

По долгу службы пришлось реализовывать много больших алгоритмов на матрицах. С незапамятных времён написал свой класс вроде и к нему кучу перегруженных операторов. Разумеется, матрицы очень большие и приходится очень аккуратно работать с любыми промежуточными результатами вычислений. Первая проблема, с которой всегда сталкиваешься при написании подобной штуки — это менеджмент памяти для промежуточных вычислений. Можно городить умный огород с аллокаторами и пулом для матриц, а можно не париться и просто запертить обычные бинарные операторы, разрешив использовать только присваивающие. То есть вместо кода: Просто писал так: И всё бы ничего. С точки зрения удобства такой код меня устраивал, с точки зрения менеджмента памяти — тоже. Но отстой в том, что при исполнении первого оператора *= приходится пробегать всем элементам матрицы, потом ещё раз пробегаться по всем элементам m2 вместе с m1, а потом — ещё раз по m2 и m3. Хотя в принципе теоретически можно было бы пробегаться сразу по трём матрицам и считать в духе и сильно сэкономить на memory layout'е. Особенно заметна экономия на дешёвых операциях вроде умножения матрицы на скаляр, но каждый раз писать код с циклами — неудобно. Особенно заметно неудобноство на сложных алгоритмах, где сразу по структуре цикла и не сообразить, что в нём происходит.

Для рабочих проектов я использую eigen и понимаю, что в нём все эти проблемы в принципе решены, но во-первых, он просто гигантский. Во-вторых, мне не нравится их философия, где слишком много решается за программиста под капотом где-то внутри перегруженных операторов. Иногда он может посчитать, что надо создавать промежуточные матрицы, иногда — нет. Конечно, это всё в принципе настраивается и контролируется, но вся эта логика раздувает библиотеку до умопомрачительных размеров. Для домашних проэктов я хотел что-то предельно простое, лишённое всех вышеописанных недостатков и обмазанное всевозможным синтаксическим сахаром, чтоб пользоваться было приятно.

Так родилась моя либа со 100% overhead-free отложенными вычислениями. Например, код из предыдущего примера записывается так:

Этот код развернётся на этапе компиляции в точно такой же код, как и в случае с ручными циклами. Обратите внимание, что в моём решении явно разделены понятия матрицы и её Proxy. Матрица содержит в себе, собственно, данные и создаётся только вручную. Операции над матрицами напрямую производить нельзя. Ни при каких обстоятельствах в промежуточных вычислениях не может случайно выделиться несколько сотен мегов памяти на здоровенную матрицу. А вот под Proxy я понимаю как раз промежуточные компайлтайм объекты, которые при присвоении резолвятся, преобразуются в единственный цикл по всем элементам и уже в нём записывают результат куда надо за один проход. Моей любимой фишкой этого подхода является то, что можно брать Proxy не на всю матрицу, а на её отдельные подматрицы. Например, такой код: запишет третью колонку произведения m3 * m1 в нулевую строку матрицы res, при этом вычислены будут только те элементы, которые, собственно, используются, то есть полное произведение матриц считаться не будет. Что скажете? У меня весь код этого добра занял около 600 строчек. Очень интересный факт, который лично меня удивил. Код генерится идентичным ручному написанию циклов в том и только том случае, если функции помечены как __forceinline. То есть просто inline, например, не катит и производительность просаживается на добрых 50%. Хотя, казалось бы, уж инлайнить-то компилятор должен уметь.

Было бы очень интересно взглянуть, если кто-то ещё из местных писал подобный велосипед. Ещё было бы очень интересно, как подобное можно реализовать на других языках программирования, потому что, как мне кажется, именно подобные задачи, связанные со сложным overhead-free синтаксическим сахаром наиболее удачно ложатся на C++.

Правка: 26 окт. 2017 18:27

Suslik
> именно подобные задачи, связанные со сложным overhead-free синтаксическим
> сахаром наиболее удачно ложатся на C++.
Да, С++ в этом месте рулит. Можно круто решить все в одну строку. А бонусом идет то, что сторонний человек глядя в этот код будет думать: "Что за говно здесь творится?".
И вот он например код, где по auto мне приехала какая то неведомая фигня, а потом хрясь хрясь, и с ней происходит снова неведомая фигня: И чтобы разбраться что за говно творится под низом - нужно лезть в исходный код библиотеки, или курить тонну манулов.
Собственно из-за этого недостатка данная проблема практически никак не решается в других ЯП.

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

Move конструкторы чем не подходят?

Правка: 26 окт. 2017 21:23

Фишка с получением третьем строки произведения впечатляет.
Но если все-таки надо умножать большие матрицы целиком, то их же быстрее не в лоб а тем хитрым рекурсивным алгоритмом из blas?

Suslik

Когда то, лет семь назад, делал такое. Кода конечно не осталось. Как оно называется и реализуется наверное пол часа вспоминал - никак.
В итоге пришлось лезть в гугл. Называется expression templates. Вот накидал пример:
http://rextester.com/BKFR72853
Основная идея - CRTP и типы на каждую операцию.

PS. Если какие косяки относительно возможностей новых стандартов, то извините, я из шаблонного "мейнстрима" выбыл еще до появления C++11.

Правка: 26 окт. 2017 21:41

старая шаблонная трава  https://github.com/highwatt/lazymatrix 

kipar
> то их же быстрее не в лоб а тем хитрым рекурсивным алгоритмом из blas?
Я бы вообще лямбду потипа for_each зафигачил для матрицы, и никаких бы прокси не изобретал. Вот это: res.GetRow(0).GetTransposed() = (m3 * m1).GetColumn(3); как по мне - слишком сложно, и нужно знать дофига инфы по конкретно этой библиотеке. Как-то так имхо проще.
Опять же непонятно как оно будет инлайнится в случае с лямбдой, и нужно тестировать производительность.
Ну и справедливости ради библиотека суслика может быть полезна в некоторых случаях, например когда с матрицами надо работать очень много.

Не матрицы и не lazy, но пусть будет
http://rextester.com/OLX73513
https://godbolt.org/g/8JTAsD

Хотел предложить expression templates, но уже предложили. Ну а так подход Суслика с view норм. Это вью называется. Для строк даже в новый стандарт завезли. Вариант с лямбдами называют комбинаторным и он не торт, ибо не читабилен. С экспрессион темплейтс можно проводить предварительные оптимизации во время компиляции. Всякие преобразования или предвычисления констант, даже посчитать обратную матрицу от константы во время компиляции. Но хоть это и торт, но компиляться будет часами даже с отключенной оптимизацией, подобных вышеописанных. Ну и отлаживать не удобно. Про мув не слушай. Чувак тебя не понял. Мув даёт экономию по памяти, но не даёт леаниризации вычислений. В общем твой вариант с вью идеален. Читается на ура даже вариант с транспонированной колонкой. Потребление памяти и линеаризация - отличные. Короче никого здесь не слушай. Все правильно сделал. Мне понравилось Пости ссылку на суслиб на гикхабе.

Правка: 26 окт. 2017 23:06

return [](){};
> Не матрицы и не lazy, но пусть будет
Ухты. Круто компилятор соптимизировал. Жаль только что assert который в default он не проверяет в компайлтайме. А то сейчас можно написать auto v2 = v1.wzya; и скомпилировать, а упадет оно только в рантайме.

*Lain*
> Вариант с лямбдами называют комбинаторным и он не торт, ибо не читабилен
Только сейчас понял, что у меня кстати вариант с лямбдами неправильный. Ибо там перемножение матриц, строка на столбец, все дела. И такое лямбдами выйдет монстроуозно. Так что признаю, библиотека вполне годная.
Тоже ссылку на гитхаб хочу. :)

MrShoor
> Ибо там перемножение матриц, строка на столбец, все дела.
Ну вот мне и непонятно. Для одной строки ок. Но для больших матриц рекурсивный алгоритм со сложностью O(n^2.8) будет сильно быстрее влобного умножения за O(n^3), а на view этот алгоритм по-моему не ложится никак (или ложится?).

А если перемножение убрать, то остается только суммирование и умножение на скаляр, возможно транспонирование - по-моему там можно исхитриться и придумать что-то без этих expression templatов, для частного случая этих трех операций. Ну, я исключительно с позиции завистника сужу.

kipar
> а на view этот алгоритм по-моему не ложится никак (или ложится?).
вот поэтому код интересно посмотреть

MrShoor
http://rextester.com/DPGKD71973