Программирование тензорной математики / Физика / Форум / Программирование игр / GameDev.ru

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

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

Рассмотрим пару примеров того, что получилось, на примере формул из этого пейпера: http://matthias-mueller-fischer.ch/publications/strainBasedDynamics.pdf
Начнём с простого, в формуле (5) вводится матрица из трёх трёхкомпонентных векторов-столбцов.

  Tensor< Index2<3, 3>, float> positions; //Index2 — второго ранга, 3, 3 — границы индексов

заполняем тензор из массива частиц:

  for (auto index : positions.Index()) //index — комбинированный индекс, состоящий из index.i(), index.j()
  {
    positions.Get(index) = particleSystem->GetParticle(triangles[index.j()])->pos[index.i()];
  }

и всё бы ничего, код на матрице выглядел бы очень похоже, но далее рассматривается случай с треугольником в трёхмерном пространстве, у которого два образующих трёхкомпонетных ребра, то есть они образуют матрицу 3x2. Куда их записать, если стандартные матрицы обычно бывают 2x2, 3x3 или 4x4? А в тензорном виде код выглядит точно так же, за исключением объявления:

  Tensor< Index2<3, 2>, float> positions; //итерирование будет автоматически осуществляться по индексам [0..2], [0..1]

обычные матричные операции реализуются несколько длиннее, чем в матричном виде, но, что очень важно, все транспонирования задаются в явном виде и ошибки вроде умножения на матрицу слева или справа допустить гораздо труднее, потому что в тензорном виде не играет роли порядок умножения, важен лишь порядок индексов. К примеру, матричная формула вида
\(S=F^TF\), записывается в тензорном виде так:
\(S_{ij}=F_{ki}F_{kj}\) и прямым образом переводится в код, в явном виде указывая суммирование по повторяющемуся индексу:

  Tensor< Index2<2, 2>, float> strain;
  for (auto index : strain.Index())
  {
    strain.Get(index) = 0;
    for (auto sum : Index1<3>())
    {
      strain.Get(index) += deformations.Get(sum.i(), index.i()) * deformations.Get(sum.i(), index.j());
    }
  }

окей, внимательный читатель может заметить, весь этот код можно реализовать на матрицах кастомного размера вроде Matrix<2, 3, float> и через обычные перегруженные операторы это бы записывалось короче:

Matrix<2, 2, float> strain = deformations * deformations.Transposed();

но вся прелесть тензорного кода в том, что он умеет то, что не умеют матрицы. Например, в Appendix A начинаются волшебные формулы (29):
\(\nabla S_{ij}=[\nabla p_{1} \nabla p_2]S_{ij}=f_j c^T_i + f_i c^T_j\) — эта формула записана в кривом постолбцовом виде, потому что в матричном виде её сформулировать невозможно и лично я потратил добрых полчаса, пытаясь понять, откуда она взялась. А в тензорном виде формула-то на самом деле простая и прозрачная:
\([\nabla_{p_l} S_{ij}]_k=\frac{dS_{ij}}{d(p_k)_l}}=(\nabla S)_{ijkl}=Q^{-1}_{li}F_{kj}+Q^{-1}_{lj}F_{ki}\)
и так же просто переводится в код:

  Tensor< Index4<2 /*strain i*/, 2/*strain j*/, 3/*gradient*/, 2/*point*/>, float > strainGradient;
  for (auto index : strainGradient.Index())
  {
    strainGradient.Get(index) =
      defPositionsInv.Get(index.l(), index.i()) * deformations.Get(index.k(), index.j()) +
      defPositionsInv.Get(index.l(), index.j()) * deformations.Get(index.k(), index.i());
  }

Вот такие дела. Интересно, что думают другие ребята, программировавшие мотан — пробовали что-то подобное? Что в моём решении можно улучшить?

}:+()___ [Smile]
Постоялец

#1
20:50, 4 мар. 2018

Suslik
> Начнём с простого, в формуле (5) вводится матрица из трёх трёхкомпонентных векторов-столбцов.
Строго говоря, такое тензором называть нельзя, ибо тензор — это набор чисел + закон их трансформации.
Но ладно, не суть.

> комбинированный индекс, состоящий из index.i(), index.j()
Лучше уж сразу index(0), index(1) — нету ограничения на размерность, да и все эти индексы по-любому немые, т. е. назначение буквы произвольно.

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

Suslik
Модератор

#2
1:31, 5 мар. 2018

}:+()___ [Smile]
> Лучше уж сразу index(0), index(1) — нету ограничения на размерность, да и все
> эти индексы по-любому немые, т. е. назначение буквы произвольно.
хоть это и удобнее с точки зрения программинга, но я к коду всегда прилагаю фотку с формулами на листочке, и если там индексы латиницей, то проще сравнивать.

> Еще хорошо бы индексы сделать типизированными, а не просто числами, чтобы была
> автоматическая проверка корректности.
вот это — очень крутая идея, я согласен. но я не придумал, как её реализовать без существенного усложнения кода. мне кажется, если реализовать эту штуку, это решит огромное количество проблем в духе умножения на матрицу слева или справа и с транспонированием квадратных матриц.

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

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

Suslik
Модератор

#3
3:20, 5 мар. 2018

}:+()___ [Smile]
в общем, сделал компайлтаймовый type safety для индексов! пример использования. сначала объявляются все размерности:

  struct WorldName {}; //просто пустые структуры, чтоб одну к другой привести нельзя было
  struct RefName {};
  struct ParticleName {};
  struct EdgeName {};
  typedef Dim<WorldName, 3> WorldDim; //мировые координаты
  typedef Dim<RefName, 2> RefDim; //материальные координаты
  typedef Dim<ParticleName, 3> ParticleDim; //номер частицы в треугольнике
  typedef Dim<EdgeName, 2> EdgeDim; //номер ребра в треугольнике, перенесённого в начало координат (первая вершина = 0)

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

  Tensor< Index2<WorldDim, RefDim>, float > deformations; //матрица деформаций. F — якобиан частных производных мировых координат(WorldDim, трёхмерные) по материальным координатам (RefDim, двумерные)
  deformations = ..; //инициализацию опустим
  Tensor< Index2<EdgeDim, RefDim>, float > defPositionsInv; //обратная матрица материальных координат, EdgeDim — столько рёбер у треугольника, RefDim — столько у них материальных координат
  defPositionsInv = ..; //тоже опустим
  Tensor< Index4<RefDim /*strain i*/, RefDim/*strain j*/, WorldDim/*gradient j*/, EdgeDim/*ref edge k*/>, float > strainGradient; //градиент тензора напряжений
  for (auto index : strainGradient.Index())
  {
    strainGradient.Get(index) =
      defPositionsInv.Get(index.l(), index.i()) * deformations.Get(index.k(), index.j()) + //каждый индекс i(), j(), k(), l() — типизированный и проверяется на соответствие в компайлтайме!
      defPositionsInv.Get(index.l(), index.j()) * deformations.Get(index.k(), index.i()); //если перепутать местами любой индекс, оно просто не скомпилится
  }

теоретически сюда можно впилить ковариантные и контрвариантные индексы, но я сам в них постоянно путаюсь :D

если писать rendering-related преобразования, например, то можно объявить мировые и view-координаты разными индексами:

  struct Dim<WordName, 3> WorldDim; //world-координаты в мировой системе координат
  struct Dim<ViewName, 3> ViewDim; //view-координаты в системе координат камеры

тогда view матрица, которая переводит из World во View будет иметь тип

Tensor<Index2<WorldName, ViewName>, float> viewMatrix;

и помножить на неё World вектор можно будет лишь единственным способом! более того, если она инвертируется транспонированием, то ещё и из view-координат можно перевести в мировые тоже только единственным способом. хотя, если она не ортонормирована, то тут как раз бы пригодились ковариантные и контрвариантные индексы :)

}:+()___ [Smile]
Постоялец

#4
4:01, 5 мар. 2018

Suslik
> если там индексы латиницей, то проще сравнивать.
А если в формуле ik без j или, вообще, αβ? А это достаточно часто, ибо разные типы индексов обозначают разными алфавитами.
Тогда уж надо как-то так:

Tensor<float, 2, 2, 3, 2> strainGradient;
for(auto [i, j, k, l] : strainGradient.dimensions())
    strainGradient(i, j, k, l) = defPositionsInv(l, i) * deformations(k, j) + defPositionsInv(l, j) * deformations(k, i);

> в общем, сделал компайлтаймовый type safety для индексов! пример использования. сначала объявляются все размерности:
Я бы сделал как-то так:

struct SpatialDimension
{
    enum { size = 3 };
    static constexpr auto x = Index<SpatialDimension>(0);
    static constexpr auto y = Index<SpatialDimension>(1);
    static constexpr auto z = Index<SpatialDimension>(2);
};

Tensor<float, SpatialDimension, SpatialDimension> matrix;

Для шаблонов важно наличие size, а остальные члены для красивого доступа к компонентам.

Еще я подумал, что нужны симметричные/антисимметричные варианты.
Т. е. чтобы можно было объявить симметричный тензор 3х3 с 6-ю компонентами так, чтобы (1,2) и (2,1) ссылались на одну ячейку памяти.

Suslik
Модератор

#5
4:14, 5 мар. 2018

}:+()___ [Smile]
> strainGradient(i, j, k, l)
я использую доступ по комбинированному Get(index) вместо Get(index.i(), index.j(), index.k(), index.l()), потому что индекс внутри представлен просто как linearOffset, а его отдельные компоненты берутся через остатки от деления. поэтому вместо того, чтобы мучить оптимизатор постоянной сборкой взад-вперёд линейного индекса в его составляющие и обратно, я просто обращаюсь по линейному индексу, где это возможно. но согласен, через tuple'ы выглядит красивее. прям иногда хочется махнуть на оптимизатор и написать такое, потому что это прикольно.

}:+()___ [Smile]
> А если в формуле ik без j или, вообще, αβ?
я просто формулу на листочке записываю, используя те же индексы, что в программе :) но вообще да, с tuple'ами, где ты сам выбираешь имена индексов, прикольно. не знаю, как это реализовать без оверхеда. но tuple'ы для индексов выглядит здорово, да.

}:+()___ [Smile]
> Еще я подумал, что нужны симметричные/антисимметричные варианты.
> Т. е. чтобы можно было объявить симметричный тензор 3х3 с 6-ю компонентами так,
> чтобы (1,2) и (2,1) ссылались на одну ячейку памяти.
да я уж тоже думал об этом.. но блин, больно много кода нового получится.

реализовал метод из пейпера выше. удивительно, но он реально заработал после второй компиляции. после первой был баг, я в тензорное представление неправильно переводил, потому что не всё хранится в нём, но сами тензорные формулы работали правильно сразу, как только скомпилились. невиданно!

}:+()___ [Smile]
Постоялец

#6
4:33, 5 мар. 2018

Suslik
> не знаю, как это реализовать без оверхеда. но прикольно, да.
Я думаю, что в большинстве случаев оптимизатор все эти циклы развернет и все индексы заменит на константные смещения.
Но чтобы наверняка, то надо убирать циклы целиком чтобы было как-то так:

TENSOR_INDEX(SpatialDimension, i);
TENSOR_INDEX(SpatialDimension, j);
TENSOR_INDEX(SpatialDimension, k);

dst(i, j) = src1(i, k) * src2(k, j);

В таком варианте можно и ко-/контравариантность проверить, и соблюдение симметрии, и оптимизировать неплохо.
Вот только шаблонокод там будет неслабый.

Suslik
Модератор

#7
4:37, 5 мар. 2018

для сравнения с альтернативными подходами. в своё время я писал разрывный метод галёркина на произвольном функциональном базисе на неструктурированных сетках, где тензор шестого ранга тензором шестого ранга погоняет. я писал в plain-c стиле, разворачивая все суммирования руками обычными циклами, по типу:

for(int cellIndex = 0; cellIndex < 4; cellIndex++)
{
  for(int valueIndex = 0; valueIndex < valuesCount; valueIndex++)
  {
    for(int dimIndex0 = 0; dimIndex0 < dimsCount; dimIndex0++)
    {
      for(int dimIndex1 = 0; dimIndex1 < dimsCount; dimIndex1++)
      {
        cellFlux[cellIndex].values[valueIndex] += incidentCells[cellIndex].functions[dimIndex0 + dimIndex1 * dimsCount].values[valueIndex] * orientation[cellIndex].refFlux[valueIndex];
      }
    }
  }
}

как можно заметить, я старался избегать многомерных массивов вроде values[cellIndex][valueIndex][dimIndex], потому что у них не понятно, в каком порядке индексы перечисляются, вместо этого для всех вложенных массивов использовал именованные поля структур. это на самом деле зарекомендовало себя как вполне неплохой способ, но есть несколько НО:
- type safety выполняется только на уровне синтаксиса в духе наименоваий: cellSomething[cellIndex] — корректно, cellSomething[valueIndex] — вероятно, баг.
- один индекс — один цикл. иногда выглядит как шестиэтажная жесть.
- memory layout фактичеки фиксирован тем, как структуры nest'ятся в друг друга. например, штуки вида cellFlux[cellIndex].values[valuesIndex].xx хранятся легко, а вот вывернуть это в структуру, транспонировав индексы в xx.cellFlux[cellIndex].values[valueIndex] — очень криво и неудобно.
- для каждого индекса нужно придумывать своё имя в коде. для двумерных массивов получается особо тупенько, имена повторяются.

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

}:+()___ [Smile]
> TENSOR_INDEX
если начинаются макросы, то всё, я даю по тормозам. ужаснота.

Suslik
Модератор

#8
7:04, 6 мар. 2018

}:+()___ [Smile]
> Tensor<float, 2, 2, 3, 2> strainGradient;
> for(auto [i, j, k, l] : strainGradient.dimensions())
strainGradient(i, j, k, l) = defPositionsInv(l, i) * deformations(k, j) + defPositionsInv(l, j) * deformations(k, i);
попытался такой способ реализовать. надо сказать, не слишком успешно. даже объявление tuple'а таким образом:

auto [i, j, k, l] =...

не догоню, как сделать. в туплах нуб. заранее объявив, можно так:

SpaceDim i;
SpaceDim j;
RefDim k;
ParticleDim k;
std::tie(i, j, k, l) = strainGradient.dimensions();

но вот чтоб прям в auto объявить — не догоню.

UPD: мдее, эта штука называется structured binding и появится только в C++17:
http://en.cppreference.com/w/cpp/language/structured_binding

werasaimon
Постоялец

#9
20:07, 6 мар. 2018

Suslik
А от когда это , векторы и матрицы перестали быть тензорами ?

P.S: и да Мат | Программирование тензорной математики ты умножаешь тензор в ковариятном базисе на транспонированый тензор в коварятном базисе , при таком умножении ты получишь тензор их сумарного ранга S(ijkk), для того чтоб выполнялось правило Энштейна , а именно сумирования по повторяющюмуся индексу , чтоб на выходе был свёрнутый тензор S(ij), у тебя в уравнении с права должен быть тензор в контервариянтном базисе , ну то бишь индексы должны быть сверху, ну или наоборот тензор в коварятном базисе надо умножать на тензор контерварятном базисе , но результат будит уже другим , так как умножения не комотутивно ! Если мы говорим о метрических(геометрических) тензорах .

Suslik
Модератор

#10
2:29, 7 мар. 2018

короче, я сдался на попытке сделать std::tuple итератор. самое близкое, что у меня получилось на pre-17 стандарте — это так:

  almost::Tensor4< RefDim /*strain i*/, RefDim/*strain j*/, WorldDim/*gradient*/, EdgeDim/*edge*/, float > strainGradient;
  for (auto index : strainGradient.GetIndex())
  {
    std::make_tuple(i, j, k, l) = index.Unpack();
    strainGradient.Get(i, j, k, l) =
      defPositionsInv.Get(l, i) * deformations.Get(k, j) +
      defPositionsInv.Get(l, j) * deformations.Get(k, i);
  }

но все индексы надо заранее объявлять, а это раздражает. остановился на методах i(), j(), k(), l(), потому что не придумал, как лучше.

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

werasaimon
> коварятном
коварятные базисы в соседнем треде, извини.

}:+()___ [Smile]
Постоялец

#11
4:02, 7 мар. 2018

Suslik
> самое близкое, что у меня получилось на pre-17 стандарте
Я бы просто сделал несколько вложенных циклов, по типу

for(auto i : dim<Spatial>())
    for(auto j : dim<Spatial>())
        for(auto k : dim<Spatial>())

> так вот второй вариант генерит существенно более быстрый код, хотя я надеялся, что разницы не будет.
Хмм... компилятор не осиливает развернуть циклы? Тогда да, умножение будет быстрее деления.

Suslik
Модератор

#12
18:53, 23 июля 2018

короче, много читал пейперы по мультигридам. много кода прототипировал. в итоге психанул и вернулся вот к этой затее. реализовал версию на мультииндексе. потом перевёл её на вариадик-темлпейты. посмотрел, какой код из этого генерит оптимизатор. взгрустнул. в итоге пошёл по той идее, которую смаел в последнем посте завещал:

      struct Values {}; //is just used as a unique type to allow compile-time type safety
      using ValuesDim = almost::StaticDimension<Values, 3>;
      struct NullspaceVectors {};
      using NullspaceVectorsDim = almost::StaticDimension<NullspaceVectors, 7>;
      struct Basis {};
      using BasisDim = almost::DynamicDimension<Basis>;

      ValuesDim valuesDim;
      NullspaceVectorsDim nullspaceDim;
      BasisDim basisDim(1); //size specified in runtime


      auto valueToNullspace = almost::MakeStaticTensor<float, ValuesDim, NullspaceVectorsDim>(); //on stack, static array
      auto nullspaceToBasis = almost::MakeDynamicTensor<float>(nullspaceDim, basisDim); //basisDim is dynamic, allocate on heap
      auto valuesInBasis = almost::MakeDynamicTensor<float>(valuesDim, basisDim); //basisDim is dynamic, allocate on heap

      for (auto i : valuesDim)
        for (auto k : basisDim)
          valuesInBasis.Get(i, k) = (i.value == k.value) ? 1.0f : -42.0f;

      for (auto i : valuesDim)
        for (auto j : nullspaceDim)
          for (auto k : basisDim)
            valuesInBasis.Get(i, k) += valueToNullspace.Get(i, j) * nullspaceToBasis.Get(j, k);

фичи:
- один из аргументов любого типа Dimension'а — это некоторая уникальная структура, что-то вроде его названия. таким образом гарантируется, что даже одинаковые по структуре индексы, соответствующие разным названиям, присвоить нельзя на этапе компиляции.
- поддерживаются статические и динамические Dimension'ы. если у тензора все размерности статические, то его можно создать на стеке. иначе создастся на куче на основе обыкновенного std::vector'а.
- полная типобезопасность при обращении к тензору по индексам. если Get(i, j) скомпилилось, но индексы имеют разный тип, то Get(j, i) не скомпилится.
- оптимизатор счастлив, код векторизуется. на мультииндексах ему было плохо, почти ничего не векторизовалось.