Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Update III  Mathematica 10.2.0 now ships with a predefined System`Permanent function, which the PermanentCode package replaces with the compatible function PermanentCode`Permanent.

For large MachineNumber arrays (both real and complex), the new PermanentCode`Permanent is ~1000× faster than the predefined System`Permanent.

For other array types—including symbolic matrices and extended precision arrays—the results are identical and the speed is comparable.

Update II  Sample code for simulating boson-sampling experiments has been added (as an answer).

This code exploits new Mathematica capabilities relating to both empirical and smooth distributions; in particular KolmogorovSmirnovTest[__] finds use.

Update I  Multiple optimizations that were suggested by members "ssch" and Simon Woods have in aggregate yielded a ~5× code-speedup; and these optimizations now are incorporated in the example code.

Thank you both very much.

Further improvements are welcome, needless to say. In particular, for n×n matrix arguments, a further $O(n)$ speedup can be achieved (in principle) by exploiting the Gray code structure of δPermutationList. However, this would come at the cost of substantially increased code complexity and generally larger round-off error.


For research in BosonSampling (for example) it is desirable to compute matrix permanents by the fastest feasible algorithm. The appended Mathematica code uses Glynn's formula to compute the complex-valued matrix permanent. This code computes the permanent of a 20×20 matrix in ~250 ms (on a 2.93 GHz MacBook Pro laptop)

The Question Asked  Can further speed gains be achieved in numerical computation of the (complex-valued) matrix permanent?

The matrices of interest typically have dimension 10×10 to 25×25, and speed-of-execution for repeated permanent evaluations at fixed matrix-dimension is the sole figure-of-merit.

Suggestions for improvements will cheerfully be adopted!

--- code follows ---


BeginPackage["PermanentCode`"];

(* If the symbol System`Permanent exists (it was introduced circa 
   Mathematica 10.2.0) then System`Permanent is Unprotected[_],
   System`Permanent is Removed[_]; and finally, a new function 
   PermanentCode`Permanent[_] is defined.

   In general PermanentCode`Permanent[mArg] returns the same result as 
   System`Permanent[mArg].  However, numerical arguments having precision 
   MachinePrecision are evaluated by compiled C-code; for large matrices
   this C-code is 1000X (or more) faster than System`Permanent. *)

If[(* this version of Mathematica defines System`Permanent[_] *)
    "System`Permanent"//NameQ, 
    (* then remove System`Permanent[_] and issue a warning *)
    Permanent::removed = 
        "(caveat) the function `1`Permanent[_] has been removed; "<>
        "the new function `2`Permanent[_] compatibly replaces it "<>
        "(with faster evaluation).";
    Message[
        Permanent::removed,
        Permanent//Context,
        $ContextPath//First
    ];
    If[Permanent//Attributes//MemberQ[#,Protected]&,
        Unprotect[Permanent];];
    ClearAll["System`Permanent"];
    Remove["System`Permanent"];
];

(* ClearAll definitions in the present Context *)
ClearAll[Context[]<>"*"//Evaluate];

Permanent::usage = "\<\
Permanent[mArg_List?MatrixQ] is computed by Glynn's formula.  The 
algorithm requires O(m^2 2^m) operations, where m is the dimension 
of the matrix arg.

Compiled evaluation is applied solely to arguments \"mArg\" that match 
either of following patterns:

     MatrixQ[mArg,IntegerQ]
     MatrixQ[mArg,MachineNumberQ]

NEAT EXAMPLES: Integer-to-Real conversions commonly evaluate more 
quickly than \"Permanent[mArg_?MatrixQ]\", per the following idiom:

     Permanent[mArg_?Integer] := Permanent[mArg//
         SetPrecision[#,MachinePrecision]&]//Round

POSSIBLE ISSUES: For Integer arguments, the compiled C-code uses 8-byte 
integers(apparently); hence too-large integer-valued permanents elicit 
an overflow (?) error as follows: 

  CompiledFunction::cfne: Numerical error encountered; 
   proceeding with uncompiled evaluation.

NOTES 
(1) Glynn's formula is re-ordered with a view to speed-by-simplicity 
    (at negligible cost in formal efficiency); in brief the algorithm
    is implemented as a sequence of BLAS-compatible calls to built-in 
    Mathematica (BLAS) functions.
(2) At present the algorithm does not fully exploit the Gray-code 
    structure of an (internal) permutation list.

RESOURCES
URL: http://en.wikipedia.org/wiki/Computing_the_permanent#Glynn_formula
URL: http://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms
URL: http://mathematica.stackexchange.com/q/38177\>";

directPermanent::usage = "\<\
directPermanent[_] is computed (inefficiently) by \
a \"no-tricks\" combinatorical sum.\>";

Begin["`Private`"];
ClearAll[Context[]<>"*"//Evaluate];

directPermanent[\[DiamondSuit]mArg_List?SquareMatrixQ] := Module[
    {\[DiamondSuit]rowList,\[DiamondSuit]colPerms},
    \[DiamondSuit]rowList = \[DiamondSuit]mArg//Length//Range;
    \[DiamondSuit]colPerms = \[DiamondSuit]rowList//Permutations;
    Map[
        (MapThread[\[DiamondSuit]mArg[[#1,#2]]&,{\[DiamondSuit]rowList,#}]//
          Times@@#&)&,
        \[DiamondSuit]colPerms
    ]//Plus@@#&//
      (\[DiamondSuit]rowList=.;\[DiamondSuit]colPerms=.;#)&
];

(* this is Permanent's sole DownValue, i.e, Permanent 
   is defined solely as a wrap around \[DiamondSuit]Permanent *)
Permanent[\[DiamondSuit]mArg_?SquareMatrixQ] := \[DiamondSuit]Permanent[\[DiamondSuit]mArg];

(* 
-------------------------------------------
Remarks upon Precision and MachinePrecision
-------------------------------------------
The function Precision treats the symbol MachinePrecision 
in a special way: "If x is not a number, Precision[x] 
gives the minimum value of Precision for all the numbers 
that appear in x. MachinePrecision is considered smaller 
than any explicit precision." 

That is why

    {
        1.0,
        1.0//SetPrecision[#,0.5*MachinePrecision]&
    }//Precision//Print;

prints "MachinePrecision".  It follows that patterns
that match low-precision matrices have to examine the 
matrix elements individually (as below).
*)

\[DiamondSuit]Permanent[ (* low-precision evaluation wrapper *)
    \[DiamondSuit]mArg_?MatrixQ/;MemberQ[\[DiamondSuit]mArg,_?(Precision[#]<MachinePrecision&),{2}]
] := Module[{\[DiamondSuit]precision},
    \[DiamondSuit]precision = \[DiamondSuit]mArg//Precision;
    \[DiamondSuit]mArg//
      SetPrecision[#,MachinePrecision]&//
        \[DiamondSuit]Permanent//
          SetPrecision[#,\[DiamondSuit]precision]&
];

\[DiamondSuit]glynnSignList::usage = "\<\
List of Gray-code permutations, saved in-memory 
for use by Permanent[_]'s Glynn-formula.\>";

\[DiamondSuit]glynnSignListMostRecentArgument = 1;

(* ensure that the only DownValues stored for
   \[DiamondSuit]glynnSignList[\[DiamondSuit]mArg_Integer] are for
   \[DiamondSuit]mArg = 1 and \[DiamondSuit]mArg = \[DiamondSuit]glynnSignListMostRecentArgument *)

\[DiamondSuit]glynnSignList[1] := (
    If[\[DiamondSuit]glynnSignListMostRecentArgument>1,
        \[DiamondSuit]glynnSignList[
          \[DiamondSuit]glynnSignListMostRecentArgument]=.;
        \[DiamondSuit]glynnSignListMostRecentArgument = 1;
    ];
    {{1}}
);

\[DiamondSuit]glynnSignList[\[DiamondSuit]m_Integer]/;(\[DiamondSuit]m>1) := (
    \[DiamondSuit]glynnSignList[\[DiamondSuit]m] = \[DiamondSuit]glynnSignList[\[DiamondSuit]m-1]//
      (   (* to conserve memory, purge unneeded DownValues *) 
          If[(\[DiamondSuit]m-1)>1,\[DiamondSuit]glynnSignList[\[DiamondSuit]m-1]=.;];
          If[\[DiamondSuit]m<\[DiamondSuit]glynnSignListMostRecentArgument,
              \[DiamondSuit]glynnSignList[
                  \[DiamondSuit]glynnSignListMostRecentArgument]=.;];
          \[DiamondSuit]glynnSignListMostRecentArgument = \[DiamondSuit]m;
          #
      )&//( Map[({1,1}~Join~(#//Rest))&,#] ~ Join ~
        Map[({1,-1}~Join~(#//Rest))&,#//Reverse] )&
);

\[DiamondSuit]compiledGlynnProductInteger = Compile[{
        {\[DiamondSuit]d, _Integer, 1},
        {\[DiamondSuit]a, _Integer, 2}
    },
    Apply[Times,(\[DiamondSuit]d.\[DiamondSuit]a)],
    CompilationTarget -> "C", 
    RuntimeAttributes -> {Listable},
    Parallelization -> True
    (*
    (* Caveat: enable for ~2x speed, less robustness *)
    ,RuntimeOptions -> {CatchMachineIntegerOverflow ->False}
    *)
];

\[DiamondSuit]compiledGlynnProductReal = Compile[{
        {\[DiamondSuit]d, _Integer, 1},
        {\[DiamondSuit]a, _Real, 2}
    },
    Apply[Times,(\[DiamondSuit]d.\[DiamondSuit]a)],
    CompilationTarget -> "C", 
    RuntimeAttributes -> {Listable},
    Parallelization -> True
];

\[DiamondSuit]compiledGlynnProductComplex = Compile[{
        {\[DiamondSuit]d, _Integer, 1},
        {\[DiamondSuit]a, _Complex, 2}
    },
    Apply[Times,(\[DiamondSuit]d.\[DiamondSuit]a)],
    CompilationTarget -> "C",
    RuntimeAttributes -> {Listable},
    Parallelization -> True
];

(* 
----------------------------------------------
Remarks upon "RuntimeAttributes -> {Listable}" 
----------------------------------------------

For compiled functions Mathematica applies RuntimeAttribute 
"Listable" attribute differently than for rule-based functions; 
namely: "When the arguments [of a 'Listable' compiled function]
include a list with higher rank than the input specification, 
the function threads over that argument."

See: Compile/tutorial/Operation#76381003

Thus we have 

  f = Compile[{{a, _Integer, 1},{b, _Integer, 2}},
      {a,b//Flatten},
      CompilationTarget -> "C", 
      RuntimeAttributes -> {Listable},
      Parallelization -> True
  ];
  f[{{1,2}},{{3,4}}] === {{{1, 2}, {3, 4}}}

whereas in contrast, a non-compiled version of the same 
Listable function threads over *all* arguments

  SetAttributes[g,Listable]
  g[a_,b_] := {a,b};
  g[{{1,2}},{{3,4}}] === {{{1, 3}, {2, 4}}} 

The following Permanent//DownValues relies crucially upon the 
just-described "RuntimeAttributes -> {Listable}" behavior 
of compiled functions.
*)

\[DiamondSuit]Permanent[ (* purely _Integer matrices *)
    \[DiamondSuit]mArg_?(MatrixQ[#,IntegerQ]&)
] := \[DiamondSuit]compiledGlynnProductInteger[
    \[DiamondSuit]glynnSignList[\[DiamondSuit]mArg//Length],
    \[DiamondSuit]mArg
]//Total[#[[1 ;; ;; 2]]] - Total[#[[2 ;; ;; 2]]]&//
  #/2^((\[DiamondSuit]mArg//Length)-1)&;

\[DiamondSuit]Permanent[ (* purely _Real MachineNumber matrices *)
    \[DiamondSuit]mArg_?((MatrixQ[#,MachineNumberQ] && FreeQ[#,_Complex,{2}])&)
] := \[DiamondSuit]compiledGlynnProductReal[
    \[DiamondSuit]glynnSignList[\[DiamondSuit]mArg//Length],
    \[DiamondSuit]mArg
]//Total[#[[1 ;; ;; 2]]] - Total[#[[2 ;; ;; 2]]]&//
  #/2^((\[DiamondSuit]mArg//Length)-1)&;

\[DiamondSuit]Permanent[ (* by default, at least one _Complex MachineNumber *)
    \[DiamondSuit]mArg_?(MatrixQ[#,MachineNumberQ]&)
] := 
(* the following encompasses the general case of pure _Complex 
   "MachineNumberQ" matrices, and also mixed _Real and _Complex 
   "MachineNumberQ" matrices, by virtue of a "CoerceTensor" call
   in the compiled C-code *)
    \[DiamondSuit]compiledGlynnProductComplex[
        \[DiamondSuit]glynnSignList[\[DiamondSuit]mArg//Length],
        \[DiamondSuit]mArg
]//Total[#[[1 ;; ;; 2]]] - Total[#[[2 ;; ;; 2]]]&//
  #/2^((\[DiamondSuit]mArg//Length)-1)&;

\[DiamondSuit]Permanent[ (* the most general case; including symbolic
    extended-precision, and mixed-type matrices; thus 
    including (for example) matrix arguments that match 
       (MemberQ[#,_Integer,{2}] && MemberQ[#,_Real,{2}])&
    and hence match no prior \[DiamondSuit]Permanent DownValue. *) 
    \[DiamondSuit]mArg_?MatrixQ
] := Map[
    Apply[Times,#.\[DiamondSuit]mArg]&,
    \[DiamondSuit]glynnSignList[\[DiamondSuit]mArg//Length]
]//Total[#[[1 ;; ;; 2]]] - Total[#[[2 ;; ;; 2]]]&//
  #/2^((\[DiamondSuit]mArg//Length)-1)&;

End[];
EndPackage[];

Code to validate and benchmark

nPerm = 4;
Table[\[DoubleStruckCapitalC][i,j],{i,1,nPerm},{j,1,nPerm}]//
  Permanent[#]-directPermanent[#]&//
    Expand//
      If[
          #===0,
          Print["VALIDATED: ",nPerm,"\[Cross]",nPerm," symbolic permanent"];,
          Print["ERROR: ",nPerm,"\[Cross]",nPerm," symbolic permanent"];
      ]&;

nPerm = 5;
Table[\[DoubleStruckCapitalC][i,j],{i,1,nPerm},{j,1,nPerm}]//
  Permanent[#]-directPermanent[#]&//
    Expand//
      If[
          #===0,
          Print["VALIDATED: ",nPerm,"\[Cross]",nPerm," symbolic permanent"];,
          Print["ERROR: ",nPerm,"\[Cross]",nPerm," symbolic permanent"];
      ]&;

nPerm = 6;
nPerm//{#,#}&//(
            1*RandomVariate[NormalDistribution[0,1],#]+
            I*RandomVariate[NormalDistribution[0,1],#]
  )*1/Sqrt[2]&//
    {Permanent[#],directPermanent[#]}&//
      (#[[1]]-#[[2]])/Sqrt[#[[2]]\[Conjugate]*#[[2]]]&//
      If[
          Abs[#]<1000*10^(-$MachinePrecision),
          Print["VALIDATED: ",nPerm,"\[Cross]",nPerm," compiled numeric permanent"];,
          Print["ERROR: ",nPerm,"\[Cross]",nPerm," compiled numeric permanent"];
      ]&;

nPerm = 7;
nPerm//{#,#}&//(
            1*RandomVariate[NormalDistribution[0,1],#]+
            I*RandomVariate[NormalDistribution[0,1],#]
  )*1/Sqrt[2]&//
    {Permanent[#],directPermanent[#]}&//
      (#[[1]]-#[[2]])/Sqrt[#[[2]]\[Conjugate]*#[[2]]]&//
      If[
          Abs[#]<100*10^(-$MachinePrecision),
          Print["VALIDATED: ",nPerm,"\[Cross]",nPerm," compiled numeric permanent"];,
          Print["ERROR: ",nPerm,"\[Cross]",nPerm," compiled numeric permanent"];
      ]&;

Do[
nPerm//{#,#}&//(
            1*RandomVariate[NormalDistribution[0,1],#]+
            I*RandomVariate[NormalDistribution[0,1],#]
  )*1/Sqrt[2]&//(
        (* first call stores Gray-code array *)
        If[iDummy==1,Print["--------------"];];
        (Permanent[#]//AbsoluteTiming)//First//1000*#&//Round//
          Print["Benchmark: ",
          Switch[iDummy,1," first",2,"second"],
          " Permanent[ ",nPerm,"\[Cross]",nPerm," ] took ",#," ms"]&;
    )&;,{nPerm,12,20},{iDummy,1,2}];
Do[
nPerm//{#,#}&//(
            1*RandomVariate[NormalDistribution[0,1],#]+
            I*RandomVariate[NormalDistribution[0,1],#]
  )*1/Sqrt[2]&//(
        (* first call stores Gray-code array *)
        If[iDummy==1,Print["--------------"];];
        (Permanent[#]//AbsoluteTiming)//First//NumberForm[#,{3,1}]&//
          Print["Benchmark: ",
          Switch[iDummy,1," first",2,"second"],
          " Permanent[ ",nPerm,"\[Cross]",nPerm," ] took ",#," s"]&;
    )&;,{nPerm,21,25},{iDummy,1,2}];

Results of validating and benchmarking

VALIDATED: 4\[Cross]4 symbolic permanent
VALIDATED: 5\[Cross]5 symbolic permanent
VALIDATED: 6\[Cross]6 compiled numeric permanent
VALIDATED: 7\[Cross]7 compiled numeric permanent
--------------
Benchmark:  first Permanent[ 12\[Cross]12 ] took 3 ms
Benchmark: second Permanent[ 12\[Cross]12 ] took 1 ms
--------------
Benchmark:  first Permanent[ 13\[Cross]13 ] took 7 ms
Benchmark: second Permanent[ 13\[Cross]13 ] took 2 ms
--------------
Benchmark:  first Permanent[ 14\[Cross]14 ] took 20 ms
Benchmark: second Permanent[ 14\[Cross]14 ] took 4 ms
--------------
Benchmark:  first Permanent[ 15\[Cross]15 ] took 33 ms
Benchmark: second Permanent[ 15\[Cross]15 ] took 8 ms
--------------
Benchmark:  first Permanent[ 16\[Cross]16 ] took 40 ms
Benchmark: second Permanent[ 16\[Cross]16 ] took 19 ms
--------------
Benchmark:  first Permanent[ 17\[Cross]17 ] took 47 ms
Benchmark: second Permanent[ 17\[Cross]17 ] took 13 ms
--------------
Benchmark:  first Permanent[ 18\[Cross]18 ] took 97 ms
Benchmark: second Permanent[ 18\[Cross]18 ] took 25 ms
--------------
Benchmark:  first Permanent[ 19\[Cross]19 ] took 198 ms
Benchmark: second Permanent[ 19\[Cross]19 ] took 51 ms
--------------
Benchmark:  first Permanent[ 20\[Cross]20 ] took 397 ms
Benchmark: second Permanent[ 20\[Cross]20 ] took 104 ms
--------------
Benchmark:  first Permanent[ 21\[Cross]21 ] took 0.8 s
Benchmark: second Permanent[ 21\[Cross]21 ] took 0.2 s
--------------
Benchmark:  first Permanent[ 22\[Cross]22 ] took 1.6 s
Benchmark: second Permanent[ 22\[Cross]22 ] took 0.4 s
--------------
Benchmark:  first Permanent[ 23\[Cross]23 ] took 3.2 s
Benchmark: second Permanent[ 23\[Cross]23 ] took 0.8 s
--------------
Benchmark:  first Permanent[ 24\[Cross]24 ] took 6.6 s
Benchmark: second Permanent[ 24\[Cross]24 ] took 1.6 s
--------------
Benchmark:  first Permanent[ 25\[Cross]25 ] took 15.3 s
Benchmark: second Permanent[ 25\[Cross]25 ] took 3.5 s

Note that the initial evaluation is slower than subsequent evaluations, because initial evaluation creates Gray-code tables that are retained for subsequent use.

share|improve this question
    
John, if I may: I think your package is better posted as an answer to your question than as an update. Nevertheless, it is both delightful and depressing that your permanent function is just as fast, if not faster, than the new built-in function… –  Guess who it is. Aug 4 at 0:31

5 Answers 5

up vote 16 down vote accepted

Looking at CompilePrint[compiledGlynnAlgorithm] there are some CopyTensor in it which aren't really needed. There's also a few CoerceTensor in there when it might be faster to just coerce the integer matrix once at the beginning.

By slightly adjusting the function all CopyTensor and CoerceTensor go away giving a small increase in speed:

compiledGlynnAlgorithmAlt = Compile[{
    {d, _Complex, 2}, {a, _Complex, 2}}, 
   Total@Map[Apply[Times, (#.a)*#] &, d],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True];


n = 20;
rc = RandomComplex[{-I - 1, I + 1}, {n, n}];
a = compiledGlynnAlgorithmAlt[δGrayCodeList[n], rc]; // AbsoluteTiming
b = compiledGlynnAlgorithm[δGrayCodeList[n], rc]; // AbsoluteTiming
a == b
(* {0.582192, Null} *)
(* {0.690600, Null} *)
(* True *)

Some more performance can be squeezed out by caching the resulting sign of each row in δGrayCodeList[n] the result is no longer exactly the same, but the relative difference is small:

δGrayCodeListSigns[n_] := δGrayCodeListSigns[n] = Times @@@ δGrayCodeList[n]

compiledGlynnAlgorithmKnownSign = 
  Compile[{{d, _Integer, 2}, {a, _Complex, 2}, {s, _Integer, 1}},
   s.Map[ Apply[Times, (#.a)] &, d]
   , CompilationTarget -> "C"
   , RuntimeAttributes -> {Listable}];

n = 20;
rc = RandomComplex[{-I - 1, I + 1}, {n, n}];

a = compiledGlynnAlgorithmAlt[δGrayCodeList[n], rc]; // AbsoluteTiming
b = compiledGlynnAlgorithm[δGrayCodeList[n], rc]; // AbsoluteTiming
c = compiledGlynnAlgorithmKnownSign[
      δGrayCodeList[n], rc, δGrayCodeListSigns[n]
    ]; // AbsoluteTiming

Abs[c - b]/Abs[b]

(* {0.565806, Null} *)
(* {0.614640, Null} *)
(* {0.430388, Null} *)
(* 2.49266*10^-13 *)
share|improve this answer
    
I affirm your speed-up of approximately 10% &amp; will award this the "answer" if nothing better appears. Thank you "ssh"! –  John Sidles Dec 2 '13 at 19:38
    
Your answer is "accepted". Thank you, "ssch". –  John Sidles Dec 3 '13 at 2:17
    
The code now cache's the sign, with a further speed-up that for Gray code ordering, the sign-list is a strictly alternating sequence of zeros and ones. –  John Sidles Dec 3 '13 at 10:58

You might get a speed up by restricting compiledGlynnAlgorithm to work on just one row of the Gray Code list, allowing the Listable and Parallelization to come into play. I say "might" because the speed up will depend on the details of your hardware.

Redefine compiledGlynnAlgorithm like so (note that it now takes a one dimensional list for d):

compiledGlynnAlgorithm = Compile[{{d, _Integer, 1}, {a, _Complex, 2}},
  Apply[Times, (d.a) d], 
   CompilationTarget -> "C", RuntimeAttributes -> {Listable}, Parallelization -> True]

And put the Total into Permanent

Permanent[mArg_List /; (MatrixQ[mArg, NumericQ])] := 
  Total@compiledGlynnAlgorithm[δGrayCodeList[mArg // Length], mArg] // 
    #/2^((mArg // Length) - 1) &;

a bit more speed

As ssch suggested, a little more performance can be squeezed out by exploiting the fact that the product of a given row of the Gray Code list is either 1 or -1. Furthermore, these occur alternately. So we can redefine compiledGlynnAlgorithm to remove the multiplication by d:

compiledGlynnAlgorithm = Compile[{{d, _Integer, 1}, {a, _Complex, 2}},
  Apply[Times, (d.a)],
  CompilationTarget -> "C", RuntimeAttributes -> {Listable}, Parallelization -> True]

and modify Permanent to Total the odd and even rows of the result separately:

Permanent[mArg_List /; (MatrixQ[mArg, NumericQ])] :=
 Module[{x},
  x = compiledGlynnAlgorithm[δGrayCodeList[mArg // Length], mArg];
  (Total[x[[;; ;; 2]]] - Total[x[[2 ;; ;; 2]]]) // #/2^((mArg // Length) - 1) &]

On my machine this gives about a factor of 3.5 speed increase over the original code for a 20x20 matrix.

share|improve this answer
    
Simon Woods, your optimizations are outstanding, and I have incorporated them (with thanks) into the code. Thank you very much. –  John Sidles Dec 3 '13 at 11:42

Here is a variant adapted from this MathGroup thread:

permanentC = 
  Compile[{{m, _Real, 2}}, With[{len = Length[m]}, (-1)^len*Module[
      {s = {0.}, u = 0.},
      Do[
       s = N[IntegerDigits[n, 2, len]];
       u += (-1)^Round[Total[s]]*(Times @@ (m.s)),
       {n, 2^len - 1}];
      u]], CompilationTarget -> "C"];

I checked it on the test set below.

SeedRandom[11111];
testmats = Table[RandomInteger[1, {n, n}], {n, 8, 20, 2}];

It is slightly faster than Permanent from the original post. It is also slightly wronger, so to speak. The issue is cancellation error, and for the larger dimensions it shows up in the last few places. The culprit is the Times @@ (m.s) part. Those get large and we lose digits on cancelling. Possibly there is a way to reorder things so as to avoid this numeric pitfall, but offhand I don't see it.

share|improve this answer
    
When I adapted this compiled code to accept _Complex input matrices, the resulting execution times were ~2.4X longer than with the "ssch" approach. Hence, this approach is not preferred (AFAICT). –  John Sidles Dec 2 '13 at 22:33

Edits

Per the comments, the bugfix RandomChoice[__] $\Rightarrow$ RandomSelection[__] is implemented.

Also the simulated distribution of $|permanent|^2$ is fitted to LogNormalDistribution[mu,sigma]; textual diagnostics are expanded; entropy cost is reported; links to Mathematica Distribution tutorials are provided.


Distributions in Mathematica

Having upgraded to Mathematica 10.0.0, it turns out that the various built-in symbols associated to the Wolfram language tutorials Nonparametric Statistical Distributions, Derived Statistical Distributions, and Hypothesis Tests greatly facilitate the numerical simulation of boson-sampling experiments using the above-provided PermanentCode` package.

Mathematica's built-in hypothesis test KolgomorovSmirnovTest is particularly valuable; the Wikipedia entry Kolmogorov–Smirnov test provides a good introduction.

Permanents in quantum physics

The appended Mathematica code simulates the permanent-distribution of the output photons of the experiment describe by Lund et al "Boson sampling from a Gaussian state" (see PRL 2014 and arXiv:1305.4346).

Open questions

This example is preparatory to a planned TCS StackExchange question Can Kolmogorov-Smirnov tests collapse the polynomial hierarchy?; this question is inspired by conjectures set forth by Scott Aaronson and Alex Arkhipov in "The Computational Complexity of Linear Optics" (see ToC 2013 and arXiv:1011.3245 ).

Baseline capabilities

The code provided establishes baseline computational capabilities to the simulation side of the following Permanent Entropy challenge:

The Permanent Entropy challenge  Experimentally sample a permanent distribution at lower entropy-cost than indistinguishably simulating that distribution by a classical computation.

Example code

The following code generates, on a MacBook Pro, about $12000\ \text{Joule/Kelvin}$ of wall-plug entropy in indistinguishably simulating — as assessed by Kolmogorov-Smirnov test — $k=32$ samples of $n=20$-photon scattering into a $m=400$-mode boson-sampling apparatus, whose scattering matrix is chosen Haar-randomly.

Thanks to Mathematica's built-in Distribution-handling symbols, the bulk of the code is devoted to generating nice-looking graphics.

Graphical output

permanent distributions and simulations

Note  A figure key appears at the end of the following textual output.

Textual output

... begin Haar-random matrix construction ...
... Haar-random matrix constructed, begin sampling permanents ...
boson-sampling is 000% done; elapsed time 00:00:14; 03:48:22 remaining
boson-sampling is 001% done; elapsed time 00:02:14; 03:40:45 remaining
    etc;
boson-sampling is 099% done; elapsed time 03:56:31; 00:02:27 remaining
boson-sampling is 100% done; elapsed time 03:58:58; 00:00:00 remaining
... permanent sampling done, begin computing distributions ...
... distributions done, begin CDFs ...
... CDFs done, begin k-sample simulations ...
... k-sample simulations done, plot the results (be patient) ...
... finished plotting

simulation entropy cost = (11955*Joule)/Kelvin

permanent distribution fitted to:
LogNormalDistribution[0.03810299650724991, 1.5231961079779779]

Alice's raw boson-sampling modes-detected data
{4, 12, 23, ...},{22, 24, 49, ...}
{18, 26, 35, ...},{6, 17, 59, ...}
{14, 29, 33, ...},{7, 32, 34, ...}
{4, 8, 21, ...},{44, 45, 98, ...}
{11, 31, 55, ...},{8, 44, 101, ...}
...

Bob's raw boson-sampling modes-detected data
{4, 10, 17, ...},{28, 29, 34, ...}
{74, 96, 115, ...},{11, 16, 46, ...}
{18, 21, 47, ...},{2, 16, 28, ...}
{29, 66, 77, ...},{5, 8, 20, ...}
{7, 13, 20, ...},{34, 49, 83, ...}
...

Kolmogorov-Smirnov (KS) tests ...
one-sample "A"    KS test: p = 0.348142
one-sample "B"    KS test: p = 0.824036
two-sample "AvsB" KS test: p = 0.835661

    Figure key:

       solid black: InverseCDF[EmpiricalDistribution] 
                    (smoothed) of |permanent|^2
        red dotted: InverseCDF[LogNormalDistribution] 
                    (fitted) of |permanent|^2
       gray dashed: InverseCDF[EmpiricalDistribution]
                    (smoothed) of |determinant|^2    

Mathematica code

As usual, making the graphics look nice takes the most code:

Needs["PermanentCode`"]; (* provided in question asked *)

(* ------------------------------------------------- *)
(* --- set boson-sampling simulation parameters --- *)
(* ------------------------------------------------- *)

nPhoton = 20;         (* number of photons detected *)
nSampleMax = 10^5;    (* upper bound to matrix samples; 
                         nSampleMax >= 10^5 is recommended *)
tSampleMax = 6*3600;  (* time-used upper bound in seconds *)
kSample = 32;         (* number of Kolmogorov-Smirnov samples *)

nominalPower = 250 Watt; (* nominal processor power at full load *)
nominalTemperature = 300 Kelvin; (* nominal heat-sink temperature *)  


(* ---------------------------------------- *)
(* --- construct Haar-random unitaries --- *)
(* ---------------------------------------- *)
"... begin Haar-random matrix construction ..."//Print;

mNode = nPhoton^2;
iSeed = 2^nPhoton;

SeedRandom[iSeed];
Umatrix = RandomVariate[NormalDistribution[],{mNode,mNode}] + 
    I * RandomVariate[NormalDistribution[],{mNode,mNode}]//
    SingularValueDecomposition[#,mNode]&//
    #[[1]].
      (DiagonalMatrix[RandomReal[{0,2*Pi},{mNode}]//Exp[I*#]&]).
      ConjugateTranspose[#[[3]]]&;

(* ------------------------------------------------- *)
(* --- set the scale of the median |permanent|^2 --- *)
(* ------------------------------------------------- *)
"... Haar-random matrix constructed, begin sampling permanents ..."//Print;

PessoanPostulate::usage = "\<\
Per the boson-sampling experiments of Lund et al. \"Boson sampling 
from a gaussian state\" (PRL 2014, see Figure 1), let $n$ be the 
number of photons detected among $m=n^2$ output modes.  Then for 
a Haar-distributed unitary scattering the median value of the 
squared permanent is (empirically) $2^(2-n^2/m)\\Gamma(n+1)/m^n$.\
\>";

(* the first two terms are heuristic; the remaining term  
   is a numerically-fitted improvement that in most cases 
   has negligible practical implications *)
PessoanPostulate = Gamma[nPhoton+1]/mNode^nPhoton * 
    2^(2-nPhoton^2/mNode) (* * 
    (* term commented-out *)
    2^(
        -0.6288239555022707 + 
        +0.0363892632845249*nPhoton + 
        +0.6797969693300729*nPhoton^3/mNode^2 +
        -0.2824925592014731*nPhoton^4/mNode^3
    ) *) ;

(* ------------------------------ *)
(* --- pretty-print utilities --- *)
(* ------------------------------ *)

padIntegerString[iArg_?NumberQ,nArg_Integer] := 
    NumberForm[iArg//Round,nArg,NumberPadding->{"0",""}]//
  ToString//Characters//{
        Take[#,1;;-nArg-1]//StringJoin//
          StringReplace[#,RegularExpression["^[0]*"] -> ""]&,
        Take[#,-nArg;;]//StringJoin
    }&//StringJoin;

toTime[arg_?NumberQ] := {arg}//
  ((#//First//{Quotient[#,60],Mod[#,60]}&)~Join~(#//Rest))&//
  ((#//First//{Quotient[#,60],Mod[#,60]}&)~Join~(#//Rest))&//
    Map[{":",padIntegerString[#,2]}&,#]&//
      Flatten//Rest//StringJoin@@#&;

(* ------------------------------------------------------ *)
(* --- sample combinatorically random output channels --- *)
(* ------------------------------------------------------ *)

permEstimatedRMS = Sqrt[PessoanPostulate//N];
amplitudeScale = permEstimatedRMS^(-1.0/nPhoton);

SeedRandom[iSeed+1];
{
    sample$Permanent,
        sample$Determinant,
    sample$Powerproduct,
        sample$Rows,
    sample$Columns
} = For[
    iSample=0;
    lastCalibration = iSample;
    startTime = zeroTime = AbsoluteTime[];,
    iSample<nSampleMax,
    iSample++,
    rowIndexList = RandomSample[Range[mNode],nPhoton]//Sort;
    colIndexList = RandomSample[Range[mNode],nPhoton]//Sort;
    Umatrix[[rowIndexList,colIndexList]]//
        (* from a superabundance of caution, rescale 
           such that the computed permament is \[ScriptCapitalO](1) *)
        (#*amplitudeScale)&//{
            #//Permanent//#*#\[Conjugate]&//Re,
            #//Det//#*#\[Conjugate]&//Re,
            #.#\[HermitianConjugate]//Diagonal//Times@@#&//Re
        }&//#*PessoanPostulate&//
          #~Join~{rowIndexList,colIndexList}&//
            Sow;
    (* provide status messages at timely intervals *)
    If[(iSample == Max[nSampleMax/1000//Round,1]) || 
        (iSample == nSampleMax-1) ||
        ((iSample > 0) &&
          (0==Mod[iSample,nSampleMax/100//Round])),
        statusString = "boson-sampling is " <>
            (100*iSample/nSampleMax//padIntegerString[#,3]&) <> 
            "% done; elapsed time " <>
            ((AbsoluteTime[]-zeroTime)//toTime) <> "; " <>
            (((AbsoluteTime[]-startTime)/(iSample-lastCalibration+1))*
                (nSampleMax-(iSample+1))//toTime) <> " remaining";
        (* uncomment to periodically pipe status to /tmp *)
        "printf '%s\\n' \"" <> statusString <> 
            "\" >/tmp/permanentStatus.txt"//Run;
        If[ (0 == Mod[iSample,nSampleMax/20//Round]) ||
            (iSample <= nSampleMax/30) || 
            ((nSampleMax-iSample) <= nSampleMax/30),
            statusString//Print;
        ];
        lastCalibration = iSample;
        startTime = AbsoluteTime[];
    ];
]//Hold//
  TimeConstrained[#//ReleaseHold,tSampleMax]&//
    Reap//Last//Last//
      (* Sort[#]& yields increasing permanent *) 
      Sort[#]&//Transpose[#]&;

samplePermanentNormed = sample$Permanent//
  Map[(1/PessoanPostulate)*#&,#]&;

sampleDeterminantNormed = sample$Determinant//
  Map[(1/PessoanPostulate)*#&,#]&;

(* ----------------------------------------------------------- *)
(* --- construct distributions both empirical and smoothed --- *)
(* ----------------------------------------------------------- *)
"... permanent sampling done, begin computing distributions ..."//Print;

(* Note: all distributions are of Log[normed |perm|^2] *)

empiricalPermanent$D = {samplePermanentNormed,samplePermanentNormed}//
  Rule[#[[1]],(#[[2]]//Log[#]&)]&//
    EmpiricalDistribution;

(* bootstrap a smooth distribution; see e.g. Wikipedia's discussion 
   URL: http://en.wikipedia.org/wiki/Bootstrapping_%28statistics%29 *)
smoothPermanent$D = empiricalPermanent$D//
  RandomVariate[#,{10*(samplePermanentNormed//Length)//Round}]&//
    SmoothKernelDistribution[#,0.1]&;

empiricalDeterminant$D = {samplePermanentNormed,sampleDeterminantNormed}//
  Rule[#[[1]],(#[[2]]//Log[#]&)]&//
    EmpiricalDistribution;

(* bootstrap a smooth distribution *)
smoothDeterminant$D = empiricalDeterminant$D//
  RandomVariate[#,{10*(sampleDeterminantNormed//Length)//Round}]&//
    SmoothKernelDistribution[#,0.05]&;

weightedPermanentData = samplePermanentNormed//
      WeightedData[#//Log[#]&,#]&;

fittedPermanent$D = weightedPermanentData//
  Module[{\[Mu]\[FilledDiamond], \[Sigma]\[FilledDiamond]},
      EstimatedDistribution[#,
          NormalDistribution[\[Mu]\[FilledDiamond], \[Sigma]\[FilledDiamond]]]
  ]&;

(* ------------------------------------------------ *)
(* --- construct inverse distributions and CDFs --- *)
(* ------------------------------------------------ *)
"... distributions done, begin CDFs ..."//Print;

permanent$CDF = smoothPermanent$D//CDF;
permanent$InverseCDF = smoothPermanent$D//InverseCDF;
fittedPermanent$InverseCDF = fittedPermanent$D//InverseCDF;
determinant$InverseCDF = smoothDeterminant$D//InverseCDF;

(* ------------------------------------------------------ *)
(* --- simulate k-sample experiments by Alice and Bob --- *)
(* ------------------------------------------------------ *)
"... CDFs done, begin k-sample simulations ..."//Print;

SeedRandom[iSeed+2];

CDFList::usage = "\<\
CDFList is a list of Rules -- intended for use with Nearest[_] --
that satisfies the relations

    sample$Permanent[[ CDFList[[i,2]] ]] \[TildeTilde] 
           InverseCDF[empiricalPermanent$D][ CDFList[[i,1]] ]

The CDFList makes it easy to simulate boson-sampling experiments.\
\>";

CDFList = sample$Permanent//
  FoldList[Plus,#]&//
    Times[#,1/(#//Last)]&//
      MapIndexed[Rule,#]&;

(* Simulate Alice's raw data *)
AlicesRawData = RandomReal[{0,1},kSample]//
  Map[(Nearest[CDFList,#]//Flatten//First)&,#]&//
    {sample$Rows[[#]],sample$Columns[[#]]}&//
      Transpose;

(* Compute Alice's processed data *)
AlicesLogNormedPermanentData = AlicesRawData//
  Map[
      (Umatrix[[ #[[1]] , #[[2]] ]]//
        (#*amplitudeScale)&//
          Permanent//#*#\[Conjugate]&//Re//Log)&,#
  ]&;

simulationEntropy = (AbsoluteTime[]-zeroTime) *
    (nominalPower/Watt)/(nominalTemperature/Kelvin)//
   Round//#*Joule/Kelvin&;

(* Simulate Bob's raw data *)
BobsRawData = RandomReal[{0,1},kSample]//
  Map[(Nearest[CDFList,#]//Flatten//First)&,#]&//
    {sample$Rows[[#]],sample$Columns[[#]]}&//
      Transpose;

(* Compute Bob's processed data *)
BobsLogNormedPermanentData = BobsRawData//
  Map[
      (Umatrix[[ #[[1]] , #[[2]] ]]//
        (#*amplitudeScale)&//
          Permanent//#*#\[Conjugate]&//Re//Log)&,#
  ]&;

(* ----------------------------- *)
(* --- plot it all up nicely --- *)
(* ----------------------------- *)
"... k-sample simulations done, plot the results (be patient) ..."//Print;

nPlotPts = 1000;

range = samplePermanentNormed//Log[#]&//Sort//
  (* don't plot outliers *)
  {#[[4;;5]]//Mean,#[[-5;;-4]]//Mean}&//
    Map[permanent$CDF,#]&;

thePermanentPointList = Range[1/2,nPlotPts]/nPlotPts//
  Select[#,(#>range[[1]])&]&//
    Select[#,(#<range[[2]])&]&//
      Map[{#,permanent$InverseCDF[#]/Log[10]}&,#//N]&;

theFittedPermanentPointList = Range[1/2,nPlotPts]/nPlotPts//
  Select[#,(#>range[[1]])&]&//
    Select[#,(#<range[[2]])&]&//
      Map[{#,fittedPermanent$InverseCDF[#]/Log[10]}&,#//N]&;

theDeterminantPointList = Range[1/2,nPlotPts]/nPlotPts//
  Select[#,(#>range[[1]])&]&//
    Select[#,(#<range[[2]])&]&//
      Map[{#,determinant$InverseCDF[#]/Log[10]}&,#//N]&;

smoothPlot = {
    thePermanentPointList,
    theFittedPermanentPointList,
    theDeterminantPointList
}//
  ListPlot[#,
    PlotJoined->True,
    PlotRange->{{0,1},{-4,2}},
    PlotStyle->{
      Directive[Black,AbsoluteThickness[1.8],Opacity[1]],
      Directive[Red,AbsoluteThickness[1.8],Dotted,Opacity[1]],
      Directive[Gray,Dashed,AbsoluteThickness[1.8],Opacity[0.6]]
    },
    AspectRatio->0.9,
    AxesOrigin->{0.5,0.0},
    Ticks->{{None,None},{None,None}},
    AxesStyle->Directive[Black,AbsoluteThickness[1.2]],
    Frame->True,
    FrameStyle->Directive[Black,AbsoluteThickness[1.2]],
    FrameTicks -> {
        {{
          Range[-4,2,1],
          {"0.0001","0.001","0.01","0.1","1","10","100"}
        }//Transpose,None},
        {{
          Range[0,1,0.2],
          {"0","0.2","0.4","0.6","0.8","1"}
        }//Transpose,None}
    },
    FrameTicksStyle->Directive[Black,AbsoluteThickness[0.6],FontSize->Medium],
    GridLines -> {
        Range[0,1,0.1],
        Outer[#1+#2&,Range[-4,1,1],Range[1,10,1]//Log[10,#]&]//
          Flatten
    },
    GridLinesStyle->Directive[Black,AbsoluteThickness[0.6],Opacity[0.35]]
  ]&;

AliceBobInverseCDFPlot = {
    AlicesLogNormedPermanentData,
    BobsLogNormedPermanentData
}//
  Map[((#//EmpiricalDistribution//
        InverseCDF[#]&)[x])&,#]&//
    Map[#/Log[10]&,#]&//
  Plot[#,{x,0,1},
      Exclusions -> None, Frame -> None, GridLines -> None, Axes -> None,
      PlotPoints->400,MaxRecursion->3,
      PlotRange->{{0,1},{-4,2}},
      PlotStyle->{
          Directive[RGBColor[0.8,0.65,0.0],AbsoluteThickness[1.2],Opacity[1]],
          Directive[RGBColor[0.0,0.65,0.8],AbsoluteThickness[1.2],Opacity[1]]
      }
  ]&;

AliceBobRegionPlot = {
    AlicesLogNormedPermanentData,
    BobsLogNormedPermanentData
}//
  Map[((#//EmpiricalDistribution//
        InverseCDF[#]&)[x])&,#]&//
  Map[#/Log[10]&,#]&//
  RegionPlot[
        {
          y<#[[1]] && y>#[[2]],
          y<#[[2]] && y>#[[1]]
        },
        {x,0,1},{y,-3,3},
        Background -> None,
        Frame -> None, GridLines -> None, Axes -> None,
        PlotRange->{{0,1},{-4,2}},
        PlotPoints->200,MaxRecursion->2,
        BoundaryStyle -> None,
        PlotStyle -> {
          Directive[RGBColor[1.0,0.75,0.0],Opacity[0.125]],
          Directive[RGBColor[0.0,0.75,1.0],Opacity[0.125]]
        }
  ]&;

"\<... finished plotting\>"//Print;

Print[""];
simulationEntropy//InputForm//
  Print["simulation entropy cost = ",#]&;

Print[""];
fittedPermanent$D//
  ReplaceAll[#,NormalDistribution->LogNormalDistribution]&//
    (Print["permanent distribution fitted to:"];
     Print[#//InputForm];)&;

Print["\nAlice's raw boson-sampling modes-detected data"];
AlicesRawData//#[[1;;5,All,1;;3]]&//
  Map[Print[#[[1]]~Join~{"..."},",",#[[2]]~Join~{"..."}]&,#]&;
Print["..."];

Print["\nBob's raw boson-sampling modes-detected data"];
BobsRawData//#[[1;;5,All,1;;3]]&//
  Map[Print[#[[1]]~Join~{"..."},",",#[[2]]~Join~{"..."}]&,#]&;
Print["..."];

{AlicesRawData,BobsRawData}//
  Intersection@@#&//Length//
    If[#>0,
        Print[""];
        Print["Warning: there were ",#," simulation collisions"];
    ]&;

Print[""];
Print["Kolmogorov-Smirnov (KS) tests ..."];
{{"one-sample \"A\"   ","one-sample \"B\"   ","two-sample \"AvsB\""},{
    KolmogorovSmirnovTest[AlicesLogNormedPermanentData,smoothPermanent$D],
        KolmogorovSmirnovTest[BobsLogNormedPermanentData,smoothPermanent$D],
    KolmogorovSmirnovTest[AlicesLogNormedPermanentData,BobsLogNormedPermanentData]
}}//Transpose//
  Map[Print[#[[1]]," KS test: p = ",#[[2]]]&,#]&;

Print[""];
"\<\
Figure key:\n
   solid black: InverseCDF[EmpiricalDistribution] 
                (smoothed) of |permanent|^2
    red dotted: InverseCDF[LogNormalDistribution] 
                (fitted) of |permanent|^2
   gray dashed: InverseCDF[EmpiricalDistribution]
                (smoothed) of |determinant|^2\
\>"//Print;

Show[smoothPlot,AliceBobRegionPlot,AliceBobInverseCDFPlot,
    PlotLabel -> Style[
        "boson-sampling Kolmogorov-Smirnov analysis \n(n="<>
            (nPhoton//ToString)<>
            " photons, k=" <>
            (kSample//ToString)<>
            " detections, m=" <>
            (mNode//ToString)<>
            " modes)",
        FontSize->Medium,
        Black
    ],
    Background -> None,
    FrameLabel->{
      Style["\<\
cumulative probability\n\
(from \>" <> (samplePermanentNormed//Length//ToString) <>
" U-matrix samples)",FontSize->Medium,Black],
      Style["inverse CDF of\n|perm|^2\[Cross]2^((n^2/m)-2)\[Cross]m^n/n!",
      FontSize->Medium,Black]
    }
]//TraditionalForm
share|improve this answer
    
Followup The bugfix now is implemented: Oct 21 at 22:06 In the above code RandomChoice[__] should be RandomSelection[__] ... –  John Sidles Oct 31 '14 at 16:34

Here's my relatively compact implementation of Glynn's formula, which incorporates the Gray code optimization:

SetAttributes[GrayCode, Listable];
GrayCode[n_Integer] := BitXor[n, BitShiftRight[n]]

permanent[mat_?MatrixQ] /; Equal @@ Dimensions[mat] :=
          Module[{b = 2^(Length[mat] - 1)},
                 PadRight[{}, b, {1, -1}].(Times @@@
                 ((2 IntegerDigits[2 b - GrayCode[Range[0, b - 1]] - 1,
                                   2] - 1).mat))/b]

Here is a compiled version of the Gray-permuted Glynn formula. The code is adapted from Knuth's method for generating tuples in Gray order:

cPermanent = Compile[{{m, _Real, 2}},
                     Module[{n = Length[m], d, f, j, p, s},
                            d = Table[1, {n}]; f = Range[0, n - 1];
                            j = 0; s = 1; p = Times @@ Total[m];
                            While[j < n - 1,
                                  f[[1]] = 0; d[[n - j]] *= -1;
                                  s *= -1; p += s (Times @@ (d.m));
                                  f[[j + 1]] = f[[j + 2]]; f[[j + 2]] = j + 1;
                                  j = f[[1]]];
                            p/2^j], "RuntimeOptions" -> "Quality"];

(Modifying the compiled routine to handle complex matrices is completely straightforward.)

At least in the limited tests I did, these two perform quite well for symbolic and numerical matrices, compared to the MathWorld routine. I do not have version 10, so I can't say whether these routines are any better than the now built-in Permanent[] function.

For reference, here is an implementation of Ryser's formula, adapted from old code by Ilan Vardi:

permanent[mat_?MatrixQ] /; Equal @@ Dimensions[mat] :=
          Module[{n = Length[mat], l},
                 l = Range[2^n - 1];
                 (1 - 2 Mod[n, 2]) PadRight[{}, 2^n - 1, {-1, 1}].(Times @@@
                 Accumulate[JacobiSymbol[-1, l]
                            mat[[IntegerExponent[l, 2] + 1]]])]

It is a bit slower than the routine based on Glynn's formula.

share|improve this answer
3  
Using the matrix: m = N @ HilbertMatrix[20], your permanent is about 73 times faster than the built-in Permanent. Your cPermanent is about 493 times faster than the built-in. With option CompilationTarget -> "C" added to cPermanent we get about 775 times faster than built-in. Finally, changing "Quality" to "Speed" makes the function about 1765 times faster than built-in. These built-in functions are getting depressing. +1 –  RunnyKine Aug 16 at 23:55
    
Thanks for the tests, @RunnyKine! What about Ryser's formula? –  Guess who it is. Aug 17 at 0:39
    
Ryser's formula is about 50 times faster than built-in. –  RunnyKine Aug 17 at 1:04
1  
@Runny, thanks again! I now have mixed feelings that it's now not very hard to beat new built-in functions. You'd think that part of what you're paying for was all the hard research that went into making sure the internal algorithms are the best available… –  Guess who it is. Aug 17 at 1:09
2  
In the case of Permanent, I was specifically asked (told) to implement it but not make a huge project out of it. The general issue of whether/when to put in a new function with possibly suboptimal performance does get raised at times. Suffice it to say that there are considerations, such as getting new functionality up and running in finite time, that sometimes lead to such judgement calls. –  Daniel Lichtblau Aug 17 at 6:53

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.