Как я могу эффективно инициализировать этот разреженный массив в Mathematica? - PullRequest
Новая
       0

Как я могу эффективно инициализировать этот разреженный массив в Mathematica?

6 голосов
/ 27 октября 2011

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

Мой код для инициализации матрицы выглядит так: 

AbsoluteTiming[S = SparseArray[{{i_, 1} -> iaa[[i]],
    {i_, j_} /; Divisible[a[[j]], aa[[i]]] -> 1.}, {2*n, n}]]

Здесь n равно 4455, iaa - список вещественных чисел, a, aa - списки целых чисел. Я получаю вывод для этой строки 

{2652.014773,SparseArray[<111742>,{8910,4455}]}

Другими словами, для инициализации этой матрицы требуется 45 минут, хотя в ней всего 111 742 ненулевых записей. Для сравнения, на самом деле решение линейной программы занимает всего 17 секунд. Что дает?

Редактировать: Кроме того, кто-нибудь может объяснить, почему это занимает столько памяти, сколько он работает? Поскольку в пользовательское время этот расчет занимает менее десяти минут ... большая часть вычислительного времени тратится на перелистывание памяти.

Mathematica по какой-то причине хранит эту матрицу как не разреженную матрицу, пока строит ее? Потому что это было бы очень глупо.

Ответы [ 2 ]

10 голосов
/ 27 октября 2011

Вы, безусловно, можете сделать намного лучше. Вот код, основанный на низкоуровневом API-интерфейсе разреженного массива здесь , который я воспроизведу, чтобы сделать код автономным: 

ClearAll[spart, getIC, getJR, getSparseData, getDefaultElement, makeSparseArray];
HoldPattern[spart[SparseArray[s___], p_]] := {s}[[p]];
getIC[s_SparseArray] := spart[s, 4][[2, 1]];
getJR[s_SparseArray] := Flatten@spart[s, 4][[2, 2]];
getSparseData[s_SparseArray] := spart[s, 4][[3]];
getDefaultElement[s_SparseArray] := spart[s, 3];
makeSparseArray[dims : {_, _}, jc : {__Integer}, ir : {__Integer}, data_List, defElem_: 0] := 
    SparseArray @@ {Automatic, dims,  defElem, {1, {jc, List /@ ir}, data}};



Clear[formSparseDivisible];
formSparseDivisible[a_, aa_, iaa_, chunkSize_: 100] :=
  Module[{getDataChunkCode, i, start, ic, jr, sparseData, dims,  dataChunk, res},
    getDataChunkCode :=
      If[# === {}, {}, SparseArray[1 - Unitize@(Mod[#, aa] & /@ #)]] &@
        If[i*chunkSize >= Length[a],
           {},
           Take[a, {i*chunkSize + 1, Min[(i + 1)*chunkSize, Length[a]]}]];  
    i = 0;
    start = getDataChunkCode;
    i++;
    ic = getIC[start];
    jr = getJR[start];
    sparseData = getSparseData[start];
    dims = Dimensions[start];        
    While[True,
      dataChunk = getDataChunkCode;
      i++;
      If[dataChunk === {}, Break[]];
      ic = Join[ic, Rest@getIC[dataChunk] + Last@ic];
      jr = Join[jr, getJR[dataChunk]];
      sparseData = Join[sparseData, getSparseData[dataChunk]];
      dims[[1]] += First[Dimensions[dataChunk]];
    ];
    res = Transpose[makeSparseArray[dims, ic, jr, sparseData]];
    res[[All, 1]] = N@iaa;
    res]

Теперь вот время: 

In[249]:= 
n = 1500;
iaa = aa = Range[2 n];
a = Range[n];
AbsoluteTiming[res = formSparseDivisible[a, aa, iaa, 100];]

Out[252]= {0.2656250, Null}

In[253]:= AbsoluteTiming[
  res1 = SparseArray[{{i_, 1} :> 
  iaa[[i]], {i_, j_} /; Divisible[a[[j]], aa[[i]]] -> 1.}, {2*n, n}];]

Out[253]= {29.1562500, Null}

Итак, у нас есть 100-кратное ускорение для этого размера массива. И, конечно же, результаты одинаковы: 

In[254]:= Normal@res1 == Normal@res
Out[254]= True

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

EDIT

Код предполагает, что списки имеют правильную длину - в частности, a должен иметь длину n, тогда как aa и iaa - 2n. Таким образом, чтобы сравнить с другими ответами, тестовый код должен быть немного изменен (только для a): 

n = 500;
iaa = RandomReal[{0, 1}, 2 n];
a = Range[ n]; aa = RandomInteger[{1, 4 n}, 2 n];


In[300]:= 
AbsoluteTiming[U=SparseArray[ReplacePart[Outer[Boole[Divisible[#1,#2]]&,
a[[1;;n]],aa],1->iaa]]\[Transpose]]
AbsoluteTiming[res = formSparseDivisible[a,aa,iaa,100]]

Out[300]= {0.8281250,SparseArray[<2838>,{1000,500}]}
Out[301]= {0.0156250,SparseArray[<2838>,{1000,500}]}

In[302]:= Normal@U==Normal@res
Out[302]= True

РЕДАКТИРОВАТЬ 2

Ваш желаемый размер матрицы устанавливается на моем не очень быстром ноутбуке (M8) примерно за 3 секунды, а также при достаточно приличном использовании памяти: 

In[323]:= 
n=5000;
iaa=RandomReal[{0,1},2 n];
a=Range[ n];aa=RandomInteger[{1,4 n},2 n];
AbsoluteTiming[res = formSparseDivisible[a,aa,iaa,200]]

Out[326]= {3.0781250,SparseArray[<36484>,{10000,5000}]}
3 голосов
/ 27 октября 2011

Ваша конструкция вызывает Divisible 2*n^2 ~= 40 000 000 раз, так что она не слишком хорошо масштабируется!В зависимости от размера целых чисел это будет составлять от одной трети до половины времени.Чтобы проверить это, вам просто нужно запустить 

AbsoluteTiming[Table[Divisible[a[[j]], aa[[i]]], {i, 2*n}, {j, n}]]

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

Редактировать: Вот некоторый тестовый код, который сравнивает время (на моей машине) для различных конструкций: 

In[1]:= n = 500;
        iaa = RandomReal[{0, 1}, 2 n];
        a = Range[2 n]; aa = RandomInteger[{1, 4 n}, 2 n]; 

In[4]:= AbsoluteTiming[S = SparseArray[{{i_, 1} :> iaa[[i]], 
                     {i_, j_} /; Divisible[a[[j]], aa[[i]]] -> 1.}, {2*n, n}]]
Out[4]= {3.423123, SparseArray[<2499>, {1000, 500}]}

In[5]:= AbsoluteTiming[T = SparseArray[ReplacePart[Table[
                     Boole[Divisible[a[[i]], aa[[j]]]], {i, 1, n}, {j, 1, 2 n}],
                     1 -> iaa]]\[Transpose]]
Out[5]= {1.524575, SparseArray[<2499>, {1000, 500}]}

In[6]:= AbsoluteTiming[U = SparseArray[ReplacePart[Outer[
                     Boole[Divisible[#1, #2]]&, a[[1 ;; n]], aa], 
                     1 -> iaa]]\[Transpose]]
Out[6]= {0.916801, SparseArray[<2499>, {1000, 500}]}

In[7]:= S == T == U
Out[7]= True

Если у вас нет другой информации о целых числах в aи aa, тогда 2*n^2 тесты с использованием Divisible будут эффективно ограничивать скорость выполнения кода.

...