3
$\begingroup$

How can I improve the speed of this code? It runs okay but is quite slow for large datasets, for example, 5 million elements.

ClearAll[filterData];
filterData[dataList_] := 
 Module[{totalEntries, removedEntries, remainingEntries, 
   determinantTest}, 
  determinantTest[{x1_, x2_, x3_, x4_, x5_, x6_, x7_, x8_, x9_, x10_, 
     x11_, x12_, x13_, x14_, x15_, x16_}] := 
   Return[Det[{{1, x1, x2, x5, x4, x8}, {x2, x1, x10, x3, x11, 
        0}, {x6, x1, x12, x14, x16, x15}, {x13, x2, x4, x6, x5, 
        x6}, {x2, x1, x3, x14, x10, x15}, {x13, x5, x8, x6, x2, 
        x6}}] == 0];
  totalEntries = Length[dataList];
  removedEntries = DeleteCases[dataList, _?determinantTest];
  remainingEntries = Complement[dataList, removedEntries];
  {totalEntries, removedEntries, remainingEntries}]

(*Example usage:*)
sampleData = RandomInteger[{-1, 1}, {5*10^6, 16}];
result = filterData[sampleData]
Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ The DeleteCases / Complement thing seems unnecessary. Just use Select[dataList, determinantTest]. Or if you really need the removed/remaining thing, then there's also GroupBy[dataList, determinantTest] which gets you an Association instead of a list. $\endgroup$
    – flinty
    Commented Mar 5 at 11:23
  • $\begingroup$ @flinty, note that there is not much speed difference between the two methods ... (GatherBy can actually be slower). $\endgroup$
    – Domen
    Commented Mar 5 at 11:30
  • 3
    $\begingroup$ Also the example usage is maybe not representative of your real data. It's almost certain that a random real 16x16 matrix has non-zero determinant. Maybe RandomInteger would be a better test case. $\endgroup$
    – flinty
    Commented Mar 5 at 11:37

1 Answer 1

Reset to default
1
$\begingroup$

For more performance on the CPU you may have to look into calling into an external library or using Mathematica's compilation tools.

But if you have OpenCLLink set up and a GPU installed, here's this abomination. It's very fast and could be faster (and a lot cleaner had I not unrolled the determinant). Runtime is about ~0.5 seconds on an NVIDIA 1080 Ti.

Needs["OpenCLLink`"];
source = "
__kernel void singular(__global mint* xs, __global mint* singular_mask)
{
int i = get_global_id(0);
int x1  = xs[i * 16 + 0];
int x2  = xs[i * 16 + 1];
int x3  = xs[i * 16 + 2];
int x4  = xs[i * 16 + 3];
int x5  = xs[i * 16 + 4];
int x6  = xs[i * 16 + 5];
int x7  = xs[i * 16 + 6];
int x8  = xs[i * 16 + 7];
int x9  = xs[i * 16 + 8];
int x10 = xs[i * 16 + 9];
int x11 = xs[i * 16 + 10];
int x12 = xs[i * 16 + 11];
int x13 = xs[i * 16 + 12];
int x14 = xs[i * 16 + 13];
int x15 = xs[i * 16 + 14];
int x16 = xs[i * 16 + 15];

int d = -(x11*x14*x4*x5*x6*x8)+x10*(x14*x16*x5*x6-x14*x4*x5*x6*x6+x14*
x5*x5*x6*x8+x10*x5*(-(x14*x6)+x5*x6*x6-x6*x6*x8))+x13*(x15*(x3*(-(x3*
x4*x5)+x11*x5*x5)+x10*(-(x10*x5*x5)+x16*x5*x5)+x2*((-x10+x16)*x2*x3+
x10*(x10*x5-x16*x5)+x3*(x3*x4+x10*x5-x11*x5-x16*x5)))+x10*(x10*x14*x5*
x8-x14*x16*x5*x8)+x3*(-(x11*x14*x5*x8)+x16*x3*x5*x8)+x2*(x10*(-(x10*
x14*x8)+x14*x16*x8)+x3*(x11*x14*x8-x16*x3*x8)))+x3*(x10*x4*x5*x6*x8+
x3*(-(x16*x5*x6)+x4*x5*x6*x6-x5*x5*x6*x8)+x5*(-(x11*x5*x6*x6)+x11*(
x14*x6+x6*x6*x8)))+x12*(-(x11*x14*x5*x6)+x10*x3*x5*x6+x15*(-(x3*x5*x5)
+x11*x5*x6+x2*(x5*x5*x5-x11*x6-x4*x5*x6+x2*(x3-x2*x5+x4*x6)))+x13*(
x15*(x3*x4*x5-x11*x5*x5+x2*(-(x3*x4)+x11*x5))+x11*x14*x5*x8-x10*x3*x5*
x8+x2*(-(x11*x14*x8)+x10*x3*x8))+x2*(x11*x14*x6+x14*x4*x5*x6+x3*(-(
x10*x6)-x4*x5*x6+x5*x5*x8)+x10*x5*(-(x5*x6)+x6*x8)+x5*(-(x11*x6*x8)+
x5*(x11*x6-x14*x8))+x2*(-(x14*x4*x6)+x3*x4*x6-x11*x5*x6+x11*x6*x8+x2*(
x14*x8-x3*x8)+x10*(x5*x6-x6*x8))))+x15*(x11*x4*x5*x5*x6+x3*(-(x10*x4*
x5)+x16*x4*x5+x3*x5*x5-x11*x5*x6-x4*x4*x5*x6)+x10*(x10*x5*x6-x16*x5*
x6-x5*x5*x5*x6+x4*x5*x6*x6)+x2*(-(x11*x4*x5*x5)-x16*x4*x5*x5+x10*(x5*
x5*x5-x10*x6+x16*x6+x4*(x5*(x5-x6)-x6*x6))+x11*x5*x6*(-x6-x8)+x2*(x11*
x6*x6+x2*(x2*x3+x3*(x5-x6)-x11*x6-x16*x6+x10*(-x5+x6))+x3*(-x3-x5*x5+
x4*(-x6-x8))+x11*x5*(x6+x8)+x16*x5*(x6+x8)+x10*(x4*x6-x5*x8))+x3*(x11*
x6+x5*x5*(-x5+x6)+x10*x8-x16*x8+x4*(x4*x5+x5*x6+x6*x8))))+x2*(x11*x14*
x5*x6*x6+x11*x14*x4*x5*x8+x14*x16*x4*x5*x8+x11*x14*x6*x8*x8+x3*(-(x14*
x4*x5*x6)+x3*(x16*x6-x4*x6*x6)+x5*(x11*x6*x6+x14*x5*x8)+x11*(-(x14*x6)
-x6*x6*x8)+x10*(-(x5*x6*x6)-x6*x8*x8)+x16*(-(x4*x5*x8)+x5*(x5*x6-x6*
x8)))+x10*(x4*(x14*x6*x6+x14*x5*(x6-x8))-x14*x5*x5*x8+x10*(x14*x6-x5*
x6*x6+x6*x6*x8)+x16*(-(x14*x6)+x5*(-(x5*x6)+x6*x8)))+x2*(-(x11*x14*x5*
x6)+x11*x14*(-x6*x6-x8*x8)+x16*(-(x14*x5*x6)-x14*x8*x8)+x2*(x11*x14*
x6+x14*x16*x6+x10*x14*(-x6+x8)+x3*(-(x16*x6)-x14*x8))+x10*(-(x14*x4*
x6)+x14*x5*x6+x16*(x5*x6-x6*x8)+x14*(-(x6*x8)+x8*x8))+x3*(x14*x4*x6+
x10*x6*x6+x3*x6*x8+x16*(x6*x8+x8*x8))))+x1*(x16*(x14*x4*x6-x14*x6*x8)+
x5*(x11*x6*x6*x8-x14*x6*x8*x8)+x11*(-(x14*x6*x6*x8)-x6*x6*x8*x8)+x10*(
x14*x6*x8+x6*x6*x8*x8+x5*(-(x14*x6*x6)+x5*x6*x6-2*x6*x6*x8)+x4*(-(x14*
x6)+x5*x6*x6-x6*x6*x8))+x4*(-(x14*x4*x6*x6)-x11*x5*x6*x6+x14*x6*x6*x8+
x11*(x14*x6*x6+x6*x6*x8))+x12*(-(x14*x5*x6)+x3*x5*x6+x2*(x14*x6+x14*
x5*x6+x2*(-(x14*x6)+x3*x6)+x3*(-x6-x5*x6))+x15*(x5*x6+x2*(-x6+x2*x6-
x5*x6))+x13*(x15*(x3*x5-x5*x5+x2*(-x3+x5))+x14*x5*x8-x3*x5*x8+x2*(-(
x14*x8)+x3*x8)))+x3*(x3*x5*(-1+x6)*x6+x16*(-(x4*x6)+x6*x8)+x4*(x4*x6*
x6-x6*x6*x8)+x10*(x4*(1-x6)*x6+x6*(-x8+x6*x8))+x5*(x14*x6-x5*x6*x6+x6*
(x6*x8+x8*x8)))+x13*(x2*x3*(x14*x8-x3*x8)+x10*(x14*x4*x8-x14*x8*x8)+
x16*(-(x14*x4*x8)+x14*x8*x8)+x3*(-(x14*x5*x8)+x3*x5*x8+x16*(x4*x8-x8*
x8)+x10*(-(x4*x8)+x8*x8))+x15*(x2*x3*(x3-x5)+x16*(x4*x5-x5*x8)+x10*(-(
x4*x5)+x5*x8)+x3*(-(x3*x5)+x5*x5+x10*(x4-x8)+x16*(-x4+x8))))+x15*(x3*
x5*x6*(-1-x8)+x5*x5*x6*x8+x11*x6*x6*x8+x10*(x4*x6+x5*x6*x6-x6*x8)+x16*
(-(x4*x6)+x6*x8)+x4*(-(x11*x6*x6)+x4*x6*x6-x6*x6*x8)+x2*(x2*(x4*x5+x3*
(-x4-x6)+x10*x6-x2*x6+x5*x6+x6*x6)-x11*x6*x8+x5*(-x6*x6-x5*x8)+x16*(
x4*x6-x6*x8)+x4*(x11*x6-x4*x6-x5*x6+x6*x8)+x10*(-(x4*x6)-x5*x6+x6*(-
x6+x8))+x3*(x6+x4*x6+x5*(x6+x8))))+x2*(x16*(-(x14*x4*x6)+x14*x6*x8)+
x11*(x14*x6*x8+x6*x8*x8)+x4*(x14*x4*x6+x11*x5*x6+x11*(-(x14*x6)-x6*x8)
)+x3*(-(x14*x6)-x4*x4*x6+x3*(1-x6)*x6-x6*x6*x8+x16*(x4*x6-x6*x8)+x5*(-
(x14*x6)+x5*x6-x6*x8-x8*x8))+x5*(-(x11*x6*x8)+x14*(x6*x6+x8*x8))+x2*(-
(x14*x5*x6)-x14*x6*x6+x2*(x14*x6-x3*x6)-x14*x4*x8+x10*(-(x14*x6)+x5*
x6-x6*x8)+x3*(x14*x6+x4*x8+x6*(x6+x8)))+x10*(x14*x6*(x6-x8)+x4*(x14*
x6-x5*x6+x6*x8)+x6*(x6*x8-x8*x8)+x5*(x14*x6-x5*x6+x6*(-x6+2*x8)))));

singular_mask[i] = d == 0 ? 1 : 0;
}
";
singular = 
  OpenCLFunctionLoad[source, 
   "singular", {{_Integer, 1, "Input"}, {_Integer, 1, 
     "Output"}}, {16}];
data = RandomInteger[{-1, 1}, {5*10^6, 16}];
Timing[
 OCLresult = 
   First@singular[Flatten@data, ConstantArray[0, Length[data]], 
     Length@data];
 ]

singularvecs = Pick[data, OCLresult, 1];

reshape[x_] := {{1, x[[1]], x[[2]], x[[5]], x[[4]], x[[8]]}, {x[[2]], 
    x[[1]], x[[10]], x[[3]], x[[11]], 0}, {x[[6]], x[[1]], x[[12]], 
    x[[14]], x[[16]], x[[15]]}, {x[[13]], x[[2]], x[[4]], x[[6]], 
    x[[5]], x[[6]]}, {x[[2]], x[[1]], x[[3]], x[[14]], x[[10]], 
    x[[15]]}, {x[[13]], x[[5]], x[[8]], x[[6]], x[[2]], x[[6]]}};

(** Test Case **)
(** 
Timing[CPUresult=Boole[Det[reshape[#]]]==0&/@data;]
CPUresult==OCLresult
**)
Improve this answer
$\endgroup$
1
  • $\begingroup$ OCLresult is a binary mask of the matrices that are singular (Det = 0), and singularvecs is the result of picking out the vectors using that mask that produce singular matrices. The giant expression is the determinant of your matrix (as you've written it in your question). This won't work for real valued matrices, just integers, but can be adapted by changing types. $\endgroup$
    – flinty
    Commented Mar 8 at 20:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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