Более эффективный способ вычисления этого рекуррентного отношения в Mathematica

Question

Вербея открыл довольно интересную дискуссию о производительности стиля функционального программирования в Mathematica. Его можно найти здесь: Какой самый эффективный способ построения матриц больших блоков в Mathematica? )

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

c = Table[0, {2L+1}, {2L+1}];

c[[1, 1]] = 1;
Do[c[[i, i]] = e[[i - 1]] c[[i - 1, i - 1]], {i, 2, 2 L + 1}];
Do[c[[i, 1]] = (1 - e[[i - 1]]) c[[i - 1, 1]], {i, 2, 2 L + 1}];
Do[c[[i, j]] = (1 - e[[i - 1]]) c[[i - 1, j]] + 
    e[[i - 1]] c[[i - 1, j - 1]], {i, 2, 2 L + 1}, {j, 2, i - 1}];

Где e - некоторый внешний список. Есть ли способ, которым я мог бы написать это более эффективным способом? Кажется, я не могу найти какой-либо очевидный способ использования встроенных функций для достижения этой цели более идиоматическим и эффективным способом.

Я понимаю, что могу сделать только так много, поскольку этот код O(n^2), но у меня есть серия умножений матриц (всего около 6), которые в совокупности занимают меньше времени, чем это утверждение. Все, что я могу сделать, чтобы хоть немного ускорить это, внесло бы для меня заметные изменения во время выполнения.

Обновление: В соответствии с рекомендациями acl я попытался использовать Compile для ускорения выражений. Для сравнительно небольшого L = 600 я получаю 3,81 секунды для наивного Do[...], 1,54 секунды для простого старого Compile и 0,033 секунды для Compile[..., CompilationTarget->"C"].

Для более реалистичного размера L = 1200 время составляет 16,68, 0,605 и 0,132 для Do, Compile и Compile[.., CompilationTarget->"C"] соответственно. Я могу добиться того же ускорения на 2 порядка, что и ACL, упомянутый в его посте.

Answer 1

Попробуйте Compile. Здесь я определяю 3 функции: f, как вы ее определили, fc, скомпилированные (для какого-то байт-кода) и fcc, скомпилированные в C (посмотрите документацию о том, как исследовать сгенерированный код).

Сначала заставьте mma сообщить нам, если что-то не может быть скомпилировано:

SetSystemOptions["CompileOptions"->"CompileReportExternal"->True]

тогда функции:

ClearAll[f];
f=Function[{ell,e},
    Module[{c=Table[0,{2ell+1},{2ell+1}]},
        c[[1,1]]=1;
        Do[c[[i,i]]=e[[i-1]] c[[i-1,i-1]],{i,2,2 ell+1}];
        Do[c[[i,1]]=(1-e[[i-1]]) c[[i-1,1]],{i,2,2 ell+1}];
        Do[c[[i,j]]=(1-e[[i-1]]) c[[i-1,j]]+e[[i-1]] c[[i-1,j-1]],
            {i,2,2 ell+1},{j,2,i-1}];
        c
    ]
];


ClearAll[fc];
fc=Compile[{{ell,_Integer},{e,_Integer,1}},
    Module[{c},
        c=Table[0,{2ell+1},{2ell+1}];
        c[[1,1]]=1;
        Do[c[[i,i]]=e[[i-1]] c[[i-1,i-1]],{i,2,2 ell+1}];
        Do[c[[i,1]]=(1-e[[i-1]]) c[[i-1,1]],{i,2,2 ell+1}];
        Do[c[[i,j]]=(1-e[[i-1]]) c[[i-1,j]]+e[[i-1]] c[[i-1,j-1]],
            {i,2,2 ell+1},{j,2,i-1}];
        c
    ]
];

ClearAll[fcc];
fcc=Compile[{{ell,_Integer},{e,_Integer,1}},
    Module[{c},
        c=Table[0,{2ell+1},{2ell+1}];
        c[[1,1]]=1;
        Do[c[[i,i]]=e[[i-1]] c[[i-1,i-1]],{i,2,2 ell+1}];
        Do[c[[i,1]]=(1-e[[i-1]]) c[[i-1,1]],{i,2,2 ell+1}];
        Do[c[[i,j]]=(1-e[[i-1]]) c[[i-1,j]]+e[[i-1]] c[[i-1,j-1]],
            {i,2,2 ell+1},{j,2,i-1}];
        c
    ],
    CompilationTarget->"C",
    RuntimeOptions->"Speed"
];

ошибок нет, так что все в порядке. А теперь протестируйте (это на macbook с 2,4 ГГц Core 2 Duo, работающим от батареи):

ell=400;
e=RandomInteger[{0,1},2*ell];
f[ell,e];//Timing
fc[ell,e];//Timing
fcc[ell,e];//Timing

дает

{2.60925, Null}
{0.092022, Null}
{0.022709, Null}

так что версия, скомпилированная в C, здесь на два порядка быстрее.

Если вы изменили типы и получили ошибки компиляции, спросите.

РЕДАКТИРОВАТЬ: Если e содержит реалов, попробуйте

ClearAll[fc];
fc=Compile[{{ell,_Integer},{e,_Real,1}},
    Module[{c},
        c=Table[0.,{2ell+1},{2ell+1}];
        c[[1,1]]=1;
        Do[c[[i,i]]=e[[i-1]] c[[i-1,i-1]],{i,2,2 ell+1}];
        Do[c[[i,1]]=(1-e[[i-1]]) c[[i-1,1]],{i,2,2 ell+1}];
        Do[c[[i,j]]=(1-e[[i-1]]) c[[i-1,j]]+e[[i-1]] c[[i-1,j-1]],
            {i,2,2 ell+1},{j,2,i-1}];
        c
    ]
];

ClearAll[fcc];
fcc=Compile[{{ell,_Integer},{e,_Real,1}},
    Module[{c},
        c=Table[0.,{2ell+1},{2ell+1}];
        c[[1,1]]=1;
        Do[c[[i,i]]=e[[i-1]] c[[i-1,i-1]],{i,2,2 ell+1}];
        Do[c[[i,1]]=(1-e[[i-1]]) c[[i-1,1]],{i,2,2 ell+1}];
        Do[c[[i,j]]=(1-e[[i-1]]) c[[i-1,j]]+e[[i-1]] c[[i-1,j-1]],
            {i,2,2 ell+1},{j,2,i-1}];
        c
    ],
    CompilationTarget\[Rule]"C",
    RuntimeOptions\[Rule]"Speed"
];

вместо.

О том, как это работает, можно почувствовать, сказав

Needs["CompiledFunctionTools`"]
CompilePrint[fc]

и получение

        2 arguments
        18 Integer registers
        6 Real registers
        3 Tensor registers
        Underflow checking off
        Overflow checking off
        Integer overflow checking on
        RuntimeAttributes -> {}

        I0 = A1
        T(R1)0 = A2
        I12 = 0
        I1 = 2
        I3 = 1
        I14 = -1
        R0 = 0.
        Result = T(R2)2

1   I11 = I1 * I0
2   I11 = I11 + I3
3   I7 = I1 * I0
4   I7 = I7 + I3
5   I17 = I12
6   T(R2)2 = Table[ I11, I7]
7   I15 = I12
8   goto 13
9   I16 = I12
10  goto 12
11  Element[ T(R2)2, I17] = R0
12  if[ ++ I16 < I7] goto 11
13  if[ ++ I15 < I11] goto 9
14  R1 = I3
15  Part[ T(R2)2, I3, I3] = R1
16  I4 = I1 * I0
17  I4 = I4 + I3
18  I5 = I3
19  goto 26
20  I10 = I5 + I14
21  I8 = I10
22  R4 = Part[ T(R1)0, I8]
23  R5 = Part[ T(R2)2, I8, I8]
24  R4 = R4 * R5
25  Part[ T(R2)2, I5, I5] = R4
26  if[ ++ I5 < I4] goto 20
27  I4 = I1 * I0
28  I4 = I4 + I3
29  I5 = I3
30  goto 40
31  I10 = I5 + I14
32  I13 = I10
33  R4 = Part[ T(R1)0, I13]
34  R5 = - R4
35  R4 = I3
36  R4 = R4 + R5
37  R5 = Part[ T(R2)2, I13, I3]
38  R4 = R4 * R5
39  Part[ T(R2)2, I5, I3] = R4
40  if[ ++ I5 < I4] goto 31
41  I4 = I1 * I0
42  I4 = I4 + I3
43  I5 = I3
44  goto 63
45  I6 = I5 + I14
46  I17 = I3
47  goto 62
48  I16 = I5 + I14
49  I9 = I16
50  R4 = Part[ T(R1)0, I9]
51  R2 = R4
52  R4 = - R2
53  R5 = I3
54  R5 = R5 + R4
55  R4 = Part[ T(R2)2, I9, I17]
56  R5 = R5 * R4
57  I16 = I17 + I14
58  R4 = Part[ T(R2)2, I9, I16]
59  R3 = R2 * R4
60  R5 = R5 + R3
61  Part[ T(R2)2, I5, I17] = R5
62  if[ ++ I17 < I6] goto 48
63  if[ ++ I5 < I4] goto 45
64  Return

с именами регистров, указывающими их тип и т. Д. Вы также можете посмотреть на сгенерированный код C, если используете опцию "C" (но это немного сложнее для чтения).

Answer 2

Учитывая, что e уже определено (как упомянуто в вашем комментарии), вы можете получить первые два шага, как это:

firsttwosteps = DiagonalMatrix[FoldList[#1 * #2&, 1, e]] + 
ArrayFlatten[{{ {{0}} , {Table[0,{2L}]} }, 
{Transpose[{Rest@ FoldList[(1-#2)#1 &,1,e]}], Table[0,{2L},{2L}]}}]

То есть вы устанавливаете диагональ и первый столбец в соответствии с требованиями в двух матрицах и добавляете их. Это работает, потому что единственная зависимость второго шага от первого шага - это элемент в позиции {1,1}.

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

firsttwostepsT = DiagonalMatrix[FoldList[#1 * #2&, 1, e]] + 
ArrayFlatten[{{ {{0}} ,{Rest@ FoldList[(1-#2)#1 &,1,e]}  }, 
{ Table[0,{2L}, {1}], Table[0,{2L},{2L}]}}]

А потом:

Transpose @ FoldList[With[{cim1 = #1}, 
    Most@FoldList[(1 - #2[[1]]) #1 + #2[[1]] #2[[2]] &, 0., 
    Transpose[{Join[e, {0}], cim1}]]] &, First@firsttwostepsT, Join[e, {0}]]

В проверке, которую я сделал, это сохранило диагональ firsttwostepsT и выдает верхнюю треугольную матрицу до ее транспонирования.

Причина, по которой ваш код работает медленно, состоит в том, что вы повторяете определение одной и той же матрицы в целом, определяя отдельные части. Как правило, в Mathematica вам почти никогда не нужны циклы Do и конструкции типа ReplacePart. Почти всегда есть лучший способ.