[mg105035] Permanent Computation Efficiency

classic Classic list List threaded Threaded
10 messages Options
Reply | Threaded
Open this post in threaded view
|

[mg105035] Permanent Computation Efficiency

sunt05
Hi All,
Recently I've been dealing with one problem of permanent computation
and got puzzled about two algorithms on it.

I've find a defined function from Mathworld(http://
mathworld.wolfram.com/Permanent.html) as bellow:
Permanent[m_List] := With[{v = Array[x, Length[m]]}, Coefficient[Times
@@ (m.v), Times @@ v]]
Meanwhile, according to Laplace Expansion Theorem, we have another
function for permanent computation(expanding a permanent by the first
row):
LPermanent[m_List?(Length[#] == 1 &)] := Flatten[m] // First;
LPermanent[m_List] := Module[  {n = Length[m], v, mm},
  v = m[[All, 1]];
  mm = m[[All, 2 ;; -1]];
  (v.Table[LPermanent[Delete[mm, i]], {i, n}])
    ]

Comparison of the two functions is quite interesting but puzzling, the
former is much better than the other.
Running the two function with the same 10*10 matrix with Timing[]
appended, the former finished in 0.s while the other a much slower
779.912s!
However, according to the computation complexity analysis, the latter
is expected as a better one.

Is there something I didn't understand about?
Looking forward to your reply!
Thanks a lot!

Reply | Threaded
Open this post in threaded view
|

[mg105049] Re: [mg105035] Permanent Computation Efficiency

Leonid Shifrin
Hi,

The key observation (for me anyway) was that while Times
@@ (m.v) computes a symbolic product of Length[m] terms, each of which is a
sum of  Length[m]  terms, it does not expand this product. Were it expanded
(you may try wrapping Expand around  Times @@ (m.v) , but beware of the huge
memory consumption that this will cause), it would have been a huge number
of terms, and your expectation would probably be confirmed (hard to predict
though since symbolic addition/multiplication doesn't  really do much and
thus doesn't take much time, but OTOH it takes lots of memory to keep the
resulting expanded expression. Seems like this is still faster than your
second solution, but more memory - hungry). Now, it appears that Coefficient
is a pretty smart function which can deduce the coefficient of some power(s)
of the variable(s) without expanding the product, but rather probably using
some combinatorial arguments.

Regards,
Leonid



On Wed, Nov 18, 2009 at 4:00 AM, Sunt <[hidden email]> wrote:

> Hi All,
> Recently I've been dealing with one problem of permanent computation
> and got puzzled about two algorithms on it.
>
> I've find a defined function from Mathworld(http://
> mathworld.wolfram.com/Permanent.html) as bellow:
> Permanent[m_List] := With[{v = Array[x, Length[m]]}, Coefficient[Times
> @@ (m.v), Times @@ v]]
> Meanwhile, according to Laplace Expansion Theorem, we have another
> function for permanent computation(expanding a permanent by the first
> row):
> LPermanent[m_List?(Length[#] == 1 &)] := Flatten[m] // First;
> LPermanent[m_List] := Module[  {n = Length[m], v, mm},
>  v = m[[All, 1]];
>  mm = m[[All, 2 ;; -1]];
>  (v.Table[LPermanent[Delete[mm, i]], {i, n}])
>    ]
>
> Comparison of the two functions is quite interesting but puzzling, the
> former is much better than the other.
> Running the two function with the same 10*10 matrix with Timing[]
> appended, the former finished in 0.s while the other a much slower
> 779.912s!
> However, according to the computation complexity analysis, the latter
> is expected as a better one.
>
> Is there something I didn't understand about?
> Looking forward to your reply!
> Thanks a lot!
>
>

Reply | Threaded
Open this post in threaded view
|

[mg105063] Re: [mg105035] Permanent Computation Efficiency

Daniel Lichtblau
In reply to this post by sunt05
Sunt wrote:

> Hi All,
> Recently I've been dealing with one problem of permanent computation
> and got puzzled about two algorithms on it.
>
> I've find a defined function from Mathworld(http://
> mathworld.wolfram.com/Permanent.html) as bellow:
> Permanent[m_List] := With[{v = Array[x, Length[m]]}, Coefficient[Times
> @@ (m.v), Times @@ v]]
> Meanwhile, according to Laplace Expansion Theorem, we have another
> function for permanent computation(expanding a permanent by the first
> row):
> LPermanent[m_List?(Length[#] == 1 &)] := Flatten[m] // First;
> LPermanent[m_List] := Module[  {n = Length[m], v, mm},
>   v = m[[All, 1]];
>   mm = m[[All, 2 ;; -1]];
>   (v.Table[LPermanent[Delete[mm, i]], {i, n}])
>     ]
>
> Comparison of the two functions is quite interesting but puzzling, the
> former is much better than the other.
> Running the two function with the same 10*10 matrix with Timing[]
> appended, the former finished in 0.s while the other a much slower
> 779.912s!
> However, according to the computation complexity analysis, the latter
> is expected as a better one.
>
> Is there something I didn't understand about?
> Looking forward to your reply!
> Thanks a lot!

We had a similar discussion in house around ten years ago. (I found it
in an old email box). First, your second approach will be slower unless
you memo-ize intermediate results. Can be done as below.

permanent2[m_] /; Length[m]==1 := m[[1,1]]
permanent2[m_] := permanent2[m] = With[{mp=Drop[m,None,1]},
        Apply[Plus, Table[m[[j,1]]*permanent2[Drop[mp,{j}]], {j,Length[m]}]]]

This of course means you run a risk of memory overflow if dimensions get
too large. But if you do nothing to store intermediate results, this
method is about as slow as possible.

I will also mention that the relative speeds will depend on the type of
matrix in question, e.g. symbolic vs. numeric.

Here is a dense symbolic example, with all entries distinct and
algebraically independent.

mat[n_] := Array[m, {n,n}];

In[5]:= Timing[p9a = permanent[mat[9]];]

Out[5]= {22.7825, Null}

In[6]:= Timing[p9b = permanent2[mat[9]];]

Out[6]= {0.240964, Null}

In[10]:= Expand[p9a]===Expand[p9b]

Out[10]= True

Might get different results with numeric matrices. Or sparse matrices.

Daniel Lichtblau
Wolfram Research






Reply | Threaded
Open this post in threaded view
|

[mg105067] Re: [mg105035] Permanent Computation Efficiency

孙挺
In reply to this post by sunt05
Daniel, thank you for your suggestion!

However, I still didn't quite know the mechanism of the function Coefficient[].
And how to determine the computation complexities of the two permanent[]
functions by big O notation?
In my opinion, the complexity of permanent2[] is O(n!). Is it right?

Thanks a lot!

Sunt

On Thu, Nov 19, 2009 at 12:59 AM, Daniel Lichtblau <[hidden email]> wrote:

> Sunt wrote:
>
>> Hi All,
>> Recently I've been dealing with one problem of permanent computation
>> and got puzzled about two algorithms on it.
>>
>> I've find a defined function from Mathworld(http://
>> mathworld.wolfram.com/Permanent.html) as bellow:
>> Permanent[m_List] := With[{v = Array[x, Length[m]]}, Coefficient[Times
>> @@ (m.v), Times @@ v]]
>> Meanwhile, according to Laplace Expansion Theorem, we have another
>> function for permanent computation(expanding a permanent by the first
>> row):
>> LPermanent[m_List?(Length[#] == 1 &)] := Flatten[m] // First;
>> LPermanent[m_List] := Module[  {n = Length[m], v, mm},
>>  v = m[[All, 1]];
>>  mm = m[[All, 2 ;; -1]];
>>  (v.Table[LPermanent[Delete[mm, i]], {i, n}])
>>    ]
>>
>> Comparison of the two functions is quite interesting but puzzling, the
>> former is much better than the other.
>> Running the two function with the same 10*10 matrix with Timing[]
>> appended, the former finished in 0.s while the other a much slower
>> 779.912s!
>> However, according to the computation complexity analysis, the latter
>> is expected as a better one.
>>
>> Is there something I didn't understand about?
>> Looking forward to your reply!
>> Thanks a lot!
>>
>
> We had a similar discussion in house around ten years ago. (I found it in
> an old email box). First, your second approach will be slower unless you
> memo-ize intermediate results. Can be done as below.
>
> permanent2[m_] /; Length[m]==1 := m[[1,1]]
> permanent2[m_] := permanent2[m] = With[{mp=Drop[m,None,1]},
>        Apply[Plus, Table[m[[j,1]]*permanent2[Drop[mp,{j}]],
> {j,Length[m]}]]]
>
> This of course means you run a risk of memory overflow if dimensions get
> too large. But if you do nothing to store intermediate results, this method
> is about as slow as possible.
>
> I will also mention that the relative speeds will depend on the type of
> matrix in question, e.g. symbolic vs. numeric.
>
> Here is a dense symbolic example, with all entries distinct and
> algebraically independent.
>
> mat[n_] := Array[m, {n,n}];
>
> In[5]:= Timing[p9a = permanent[mat[9]];]
>
> Out[5]= {22.7825, Null}
>
> In[6]:= Timing[p9b = permanent2[mat[9]];]
>
> Out[6]= {0.240964, Null}
>
> In[10]:= Expand[p9a]===Expand[p9b]
>
> Out[10]= True
>
> Might get different results with numeric matrices. Or sparse matrices.
>
> Daniel Lichtblau
> Wolfram Research
>
>
>
>
>
>


Reply | Threaded
Open this post in threaded view
|

[mg105072] Re: [mg105067] Re: [mg105035] Permanent Computation Efficiency

Daniel Lichtblau
Sunt wrote:

> Daniel, thank you for your suggestion!
>
> However, I still didn't quite know the mechanism of the function Coefficient[].
> And how to determine the computation complexities of the two permanent[]
> functions by big O notation?
> In my opinion, the complexity of permanent2[] is O(n!). Is it right?
>
> Thanks a lot!
>
> Sunt

Coefficient is taking derivatives, without expanding the expression. Not
sure what is the complexity of this approach. It cannot be better than
exponential, if i recall permanent complexity correctly. Probably
factorial in dimension of the matrix.

The complexity of permanent2[] seems to be O(n!) for basic symbolic
matrices of the sort I used. I base that on some simple timings using

mat[n_] := Array[m, {n,n}];

Timing[p9b = permanent2[mat[9]];]
Timing[p10b = permanent2[mat[10]];]
Timing[p11b = permanent2[mat[11]];]

Since the sub-permenents share sub-sub-computations, and memo-ization is
used, I'm not sure what is the complexity in general. When I use integer
entries, the methods behave as though they might be better than O(n!),
and indeed permanent2 seems to get five times slower for every jump of 2
in n. If this reflects the actual complexity, the improvement over n! is
due to substantial use of memory.

Daniel Lichtblau
Wolfram Research


Reply | Threaded
Open this post in threaded view
|

[mg105077] Re: Permanent Computation Efficiency

sunt05
In reply to this post by sunt05
On Nov 19, 6:36 pm, Leonid Shifrin <[hidden email]> wrote:

> Hi,
>
> The key observation (for me anyway) was that while Times
> @@ (m.v) computes a symbolic product of Length[m] terms, each of which is a
> sum of  Length[m]  terms, it does not expand this product. Were it expanded
> (you may try wrapping Expand around  Times @@ (m.v) , but beware of the huge
> memory consumption that this will cause), it would have been a huge number
> of terms, and your expectation would probably be confirmed (hard to predict
> though since symbolic addition/multiplication doesn't  really do much and
> thus doesn't take much time, but OTOH it takes lots of memory to keep the
> resulting expanded expression. Seems like this is still faster than your
> second solution, but more memory - hungry). Now, it appears that Coefficient
> is a pretty smart function which can deduce the coefficient of some power (s)
> of the variable(s) without expanding the product, but rather probably using
> some combinatorial arguments.
>
> Regards,
> Leonid
>
>
>
> On Wed, Nov 18, 2009 at 4:00 AM, Sunt <[hidden email]> wrote:
> > Hi All,
> > Recently I've been dealing with one problem of permanent computation
> > and got puzzled about two algorithms on it.
>
> > I've find a defined function from Mathworld(http://
> > mathworld.wolfram.com/Permanent.html) as bellow:
> > Permanent[m_List] := With[{v = Array[x, Length[m]]}, Coefficient[Times
> > @@ (m.v), Times @@ v]]
> > Meanwhile, according to Laplace Expansion Theorem, we have another
> > function for permanent computation(expanding a permanent by the first
> > row):
> > LPermanent[m_List?(Length[#] == 1 &)] := Flatten[m] // First;
> > LPermanent[m_List] := Module[  {n = Length[m], v, mm},
> >  v = m[[All, 1]];
> >  mm = m[[All, 2 ;; -1]];
> >  (v.Table[LPermanent[Delete[mm, i]], {i, n}])
> >    ]
>
> > Comparison of the two functions is quite interesting but puzzling, the
> > former is much better than the other.
> > Running the two function with the same 10*10 matrix with Timing[]
> > appended, the former finished in 0.s while the other a much slower
> > 779.912s!
> > However, according to the computation complexity analysis, the latter
> > is expected as a better one.
>
> > Is there something I didn't understand about?
> > Looking forward to your reply!
> > Thanks a lot!

Thank you!

Could you please tell me what combinatorial knowledge is in need to
understand the function Coefficient?

Reply | Threaded
Open this post in threaded view
|

[mg105103] Re: [mg105035] Permanent Computation Efficiency

Daniel Lichtblau
In reply to this post by sunt05
Sunt wrote:

> Daniel, thank you for your suggestion!
>
> However, I still didn't quite know the mechanism of the function
> Coefficient[].
> And how to determine the computation complexities of the two permanent[]
> functions by big O notation?
> In my opinion, the complexity of permanent2[] is O(n!). Is it right?
>
> Thanks a lot!
>
> Sunt
> [...]

I was confused for some time about the complexity of permanent2. In the
numeric case it behaves like O(2^n), in the symbolic cases I tried it
seemed to be O(n!). I now think I understand this.

It is easy to check that O(2^n) sub-permanents get created when one
computes a permanent of an nxn from scratch (this means we erase all
stored DownValues after finishing a computation and before starting the
next one).

Here is a specific example. It shows the 2^n complexity, both in speed
and memory usage. That latter is not hard to reason from basic
principles, and indeed that also gives a measure of the expected
performance. Or at least a lower bound, as we'll see momentarily.

mat[n_] := Array[m, {n,n}];

Table[
   Clear[permanent2];
   permanent2[m_] /; Length[m]==1 := m[[1,1]];
   permanent2[m_] := permanent2[m] = With[{mp=Drop[m,None,1]},
     Apply[Plus, Table[m[[j,1]]*permanent2[Drop[mp,{j}]],
       {j,Length[m]}]]];
   {First[Timing[permanent2[mat[j]/.m[ii_,jj_]:>ii+jj^2]]],
     Length[DownValues[permanent2]]},
   {j,8,18}]

Out[12]= {{0.007999, 249}, {0.017998, 504}, {0.035994, 1015},
    {0.077988, 2038}, {0.167975, 4085}, {0.370943, 8180},
    {0.797878, 16371}, 1.74074, 32754}, {3.70044, 65521},
    {7.8638, 131056}, {16.7045, 262127}}

Now we'll try the purely symbolic case (of necessity, using smaller
dimensions).

Table[
   Clear[permanent2];
   permanent2[m_] /; Length[m]==1 := m[[1,1]];
   permanent2[m_] := permanent2[m] = With[{mp=Drop[m,None,1]},
     Apply[Plus, Table[m[[j,1]]*permanent2[Drop[mp,{j}]],
       {j,Length[m]}]]];
   {First[Timing[permanent2[mat[j]]]],
     Length[DownValues[permanent2]]},
   {j,7,11}]

Out[13]= {{0.009998, 122}, {0.037994, 249}, {0.249962, 504},
    {2.05569, 1015},  {21.8357, 2038}}

We are storing O(2^n) sub-permanents. But the timing is clearly O(n!).
This seeming mystery is resolved by checking the sizes (measured by
LeafCount) of the results. They are grwoing as N!, hence the time needed
to produce them must be as bad. Here is the code, barely modified, that
shows this.

Table[
   Clear[permanent2];
   permanent2[m_] /; Length[m]==1 := m[[1,1]];
   permanent2[m_] := permanent2[m] = With[{mp=Drop[m,None,1]},
     Apply[Plus, Table[m[[j,1]]*permanent2[Drop[mp,{j}]],
       {j,Length[m]}]]];
   {First[Timing[permanent2[mat[j]]]],
     LeafCount[permanent2[mat[j]]],
     Length[DownValues[permanent2]]},
   {j,7,11}]
 

Out[14]= {{0.011998, 53376, 122}, {0.037994, 427041, 249},
     {0.240964, 3843406, 504}, {2.08568, 38434101, 1015},
     {21.9807, 422775156, 2038}}

The upshot is that expression swell is making this purely symbolic case
show the n! behavior. In contrast, the integer case when done with
cached sub-permanents, is "merely" exponential in behavior.

Daniel Lichtblau
Wolfram Research

Reply | Threaded
Open this post in threaded view
|

[mg105369] Re: [mg105077] Re: Permanent Computation Efficiency

孙挺
In reply to this post by sunt05
Thanks all same!
The community is so nice:)

On Sun, Nov 29, 2009 at 2:39 AM, Leonid Shifrin <[hidden email]> wrote:

> Hi,
>
> Well, I wouldn't explain it better than Daniel did, and I know far less
> about it than he does. Since Coefficient takes derivatives, it does not need
> to expand the product, but can use the product formula for derivatives -
> this is I guess the reason for its efficiency in this case.
>
> Regards,
> Leonid
>
>
>
>> Could you please tell me what combinatorial knowledge is in need to
>> understand the function Coefficient?
>>
>>
>
Reply | Threaded
Open this post in threaded view
|

[mg105436] Re: Permanent Computation Efficiency

g.resta@iit.cnr.it
In reply to this post by sunt05
The best algorithm known for the computation of the permanent
in the general case is the Ryser algorithm and has a cost of
O(n 2^n) operations or O(n^2 2^n) depending on the implementation.

http://mathworld.wolfram.com/RyserFormula.html

It is quite simple to implement and it is much faster than the
recursive
algorithm you are discussing, unless the matrix is so sparse
that there are relative few submatrices whose permanent has to
be computed.

Several years ago I used it for computing the permanent
of 0-1 matrices, but I implemented it in C, where some binary tricks
can be
used.

On current hardware, the computation of the permanent
of a random NxN 0-1 matrix (with about one half of entries equal to 1)
takes less than 0.01 sec. for N=16 and about 1 second for N=24.
I do not know how efficient can be Mathematica in comparison.

giovanni resta
--
http://anagram.it : anagrams, alphametics, arithmogriphs,...

Reply | Threaded
Open this post in threaded view
|

[mg105458] Re: [mg105436] Re: Permanent Computation Efficiency

Daniel Lichtblau
[hidden email] wrote:

> The best algorithm known for the computation of the permanent
> in the general case is the Ryser algorithm and has a cost of
> O(n 2^n) operations or O(n^2 2^n) depending on the implementation.
>
> http://mathworld.wolfram.com/RyserFormula.html
>
> It is quite simple to implement and it is much faster than the
> recursive
> algorithm you are discussing, unless the matrix is so sparse
> that there are relative few submatrices whose permanent has to
> be computed.
>
> Several years ago I used it for computing the permanent
> of 0-1 matrices, but I implemented it in C, where some binary tricks
> can be
> used.
>
> On current hardware, the computation of the permanent
> of a random NxN 0-1 matrix (with about one half of entries equal to 1)
> takes less than 0.01 sec. for N=16 and about 1 second for N=24.
> I do not know how efficient can be Mathematica in comparison.
>
> giovanni resta
> --
> http://anagram.it : anagrams, alphametics, arithmogriphs,...

I think the memoized version of the recursive method is also showing
O(n*2^n) complexity for the class of inputs you indicate. Below is an
experiment to this effect.

In[15]:= permanent2[m_] /; Length[m] == 1 := m[[1, 1]]
permanent2[m_] :=
  permanent2[m] =
   With[{mp = Drop[m, None, 1]},
    Apply[Plus,
     Table[m[[j, 1]]*permanent2[Drop[mp, {j}]], {j, Length[m]}]]]

We time examples from dimensions 8 through 20.

In[3]:= data =
  Table[Timing[permanent2[RandomInteger[1, {n, n}]]], {n, 8, 20, 2}]
Out[3]= {{0.004999, 12}, {0.029996, 2062}, {0.12898, 8746},
   {0.609907, 1823066}, {2.77258, 194056418},
   {13.9019, 6047256830}, {63.5183, 347782811770}}

Now fit to the logs of the timings. This may or may not be an ideal
thing to do; I did it in order to decrease potential for dominance by
the larger terms.

In[4]:= timings = data[[All, 1]];
ldata = Log[2, timings];
ld2 = Transpose[{Range[8, 20, 2], ldata}];
ff = FindFit[ld2, n + Log[2, n] + b, b, n]
Out[7]= {b -> -18.4743}

We check the resulting fit.

In[8]:= timings - Table[2^(n + Log[2, n] + b) /. ff, {n, 8, 20, 2}]
Out[8]= {-0.000624668, 0.00187766, -0.00598802, -0.0199438,\
   -0.10674, 0.944957, 5.93199}

Seems not too bad, as fits go. So the absolute speed appears to be about
two orders of magnitude slower than what you were able to do in C, but
the overall complexity looks to be as advertised. Moreover, if I let the
power of n be a parameter to be fitted, I typically get values just a
hair larger than 1.

I also coded the Ryser formula. I was able to get a version that was on
speaking terms with Compile. It is still slower than the above, though
not by too much. It also seems to be showing roughly the right
complexity. And has the added benefit of not running up the memory tab.
I will mention that I did not try the Gray code optimization, so the
expected complexity in this case is O(n^2*2^n).

permanent3[m_] := With[{len = Length[m]},
   (-1)^len*Module[{s, t},
     Sum[
      s = IntegerDigits[n, 2, len];
      t = Total[s];
      (-1)^t*
       Product[Sum[m[[i, j]], {j, Flatten[Position[s, 1]]}], {i, len}]
      , {n, 2^len - 1}]
     ]]

In[46]:= data3 =
  Table[Timing[permanent3C[RandomInteger[1, {n, n}]]], {n, 8, 18, 2}]
Out[46]= {{0.009998, 125.}, {0.057991, 4066.}, {0.293956,
   43379.}, {1.60576, 2.62758*10^6}, {8.77567, 7.19182*10^7}, {44.9152,
    5.06598*10^10}}

My tentative conclusions are (i) the complexity seems to be as claimed,
and (ii) I did not find a particularly fast way to compute this in
Mathematica.

Daniel Lichtblau
Wolfram Research