6
$\begingroup$

As we know, QuadraticFormDistribution is the distribution of z.a.z + b.z + c for a real-valued p x p symmetric positive definite matrix a, a length-p vector b, a scalar c, and a p-dimensional multivariate normal vector z

I have a QuadraticFormDistribution expression fexp, which has a very complex form, and it can be simplified to the form z.a.z + b.z + c. Now, I want to get the matrix a, the vector b, and the scalar c from fexp. I think there are many methods for this simple problem. I always ExpandAll the expression fexp and use the Coefficient to get the matrix a, the vector b and the scalar c.

The problem is: I will deal with very large numbers of these kinds of quadratic expressions, so I want to find a faster method to deal with this problem.

Take an example, 

fexp=12 - 9 z[3] - 16 z[4] - 25 z[5] - 36 z[6] - 49 z[7] - 64 z[8] - 
81 z[9] - 100 z[10] - 121 z[11] - 144 z[12] + 
k^2 (-1 + 2 z[1] + 2 z[2] + 2 z[3] + 2 z[4] + 2 z[5] + 2 z[6] + 
2 z[7] + 2 z[8] + 2 z[9] + 2 z[10] + 2 z[11] + 2 z[12]) + 
z[1] (-1 + 54 a[1, 2] z[2] + 96 Sqrt[2] a[1, 3] z[3] + 
150 Sqrt[3] a[1, 4] z[4] + 432 a[1, 5] z[5] + 
294 Sqrt[5] a[1, 6] z[6] + 384 Sqrt[6] a[1, 7] z[7] + 
486 Sqrt[7] a[1, 8] z[8] + 1200 Sqrt[2] a[1, 9] z[9] + 
2178 a[1, 10] z[10] + 864 Sqrt[10] a[1, 11] z[11] + 
1014 Sqrt[11] a[1, 12] z[12]) + 
k (1 + 4 a[1, 1] z[1]^2 + 3 z[3] + 4 z[4] + 5 z[5] + 6 z[6] + 
7 z[7] + 8 z[8] + 9 z[9] + 10 z[10] + 11 z[11] + 12 z[12] + 
z[1] (1 + 18 a[1, 2] z[2] + 32 a[1, 3] z[3] + 50 a[1, 4] z[4] + 
   72 a[1, 5] z[5] + 98 a[1, 6] z[6] + 128 a[1, 7] z[7] + 
   162 a[1, 8] z[8] + 200 a[1, 9] z[9] + 242 a[1, 10] z[10] + 
   288 a[1, 11] z[11] + 338 a[1, 12] z[12]) + 
  2 (8 a[2, 2] z[2]^2 + 18 a[3, 3] z[3]^2 + 49 a[3, 4] z[3] z[4] + 
   32 a[4, 4] z[4]^2 + 64 a[3, 5] z[3] z[5] + 
   81 a[4, 5] z[4] z[5] + 50 a[5, 5] z[5]^2 + 
   81 a[3, 6] z[3] z[6] + 100 a[4, 6] z[4] z[6] + 
   121 a[5, 6] z[5] z[6] + 72 a[6, 6] z[6]^2 + 
   100 a[3, 7] z[3] z[7] + 121 a[4, 7] z[4] z[7] + 
   144 a[5, 7] z[5] z[7] + 169 a[6, 7] z[6] z[7] + 
   98 a[7, 7] z[7]^2 + 121 a[3, 8] z[3] z[8] + 
   144 a[4, 8] z[4] z[8] + 169 a[5, 8] z[5] z[8] + 
   196 a[6, 8] z[6] z[8] + 225 a[7, 8] z[7] z[8] + 
   128 a[8, 8] z[8]^2 + 144 a[3, 9] z[3] z[9] + 
   169 a[4, 9] z[4] z[9] + 196 a[5, 9] z[5] z[9] + 
   225 a[6, 9] z[6] z[9] + 256 a[7, 9] z[7] z[9] + 
   289 a[8, 9] z[8] z[9] + 162 a[9, 9] z[9]^2 + 
   169 a[3, 10] z[3] z[10] + 196 a[4, 10] z[4] z[10] + 
   225 a[5, 10] z[5] z[10] + 256 a[6, 10] z[6] z[10] + 
   289 a[7, 10] z[7] z[10] + 324 a[8, 10] z[8] z[10] + 
   361 a[9, 10] z[9] z[10] + 200 a[10, 10] z[10]^2 + 
   196 a[3, 11] z[3] z[11] + 225 a[4, 11] z[4] z[11] + 
   256 a[5, 11] z[5] z[11] + 289 a[6, 11] z[6] z[11] + 
   324 a[7, 11] z[7] z[11] + 361 a[8, 11] z[8] z[11] + 
   400 a[9, 11] z[9] z[11] + 441 a[10, 11] z[10] z[11] + 
   242 a[11, 11] z[
     11]^2 + (225 a[3, 12] z[3] + 256 a[4, 12] z[4] + 
      289 a[5, 12] z[5] + 324 a[6, 12] z[6] + 361 a[7, 12] z[7] + 
      400 a[8, 12] z[8] + 441 a[9, 12] z[9] + 
      484 a[10, 12] z[10] + 529 a[11, 12] z[11]) z[12] + 
   288 a[12, 12] z[12]^2 + 
   z[2] (1 + 25 a[2, 3] z[3] + 36 a[2, 4] z[4] + 
      49 a[2, 5] z[5] + 64 a[2, 6] z[6] + 81 a[2, 7] z[7] + 
      100 a[2, 8] z[8] + 121 a[2, 9] z[9] + 144 a[2, 10] z[10] + 
      169 a[2, 11] z[11] + 196 a[2, 12] z[12]))) + 
   2 (z[2] (-2 + 75 a[2, 3] z[3] + 108 Sqrt[2] a[2, 4] z[4] + 
   147 Sqrt[3] a[2, 5] z[5] + 384 a[2, 6] z[6] + 
   243 Sqrt[5] a[2, 7] z[7] + 300 Sqrt[6] a[2, 8] z[8] + 
   363 Sqrt[7] a[2, 9] z[9] + 864 Sqrt[2] a[2, 10] z[10] + 
   1521 a[2, 11] z[11] + 588 Sqrt[10] a[2, 12] z[12]) + 
   3 (49 a[3, 4] z[3] z[4] + 64 Sqrt[2] a[3, 5] z[3] z[5] + 
   100 Sqrt[2] a[4, 6] z[4] z[6] + 
   121 Sqrt[3] a[4, 7] z[4] z[7] + 169 a[6, 7] z[6] z[7] + 
   288 a[4, 8] z[4] z[8] + 196 Sqrt[2] a[6, 8] z[6] z[8] + 
   225 a[7, 8] z[7] z[8] + 169 Sqrt[5] a[4, 9] z[4] z[9] + 
   225 Sqrt[3] a[6, 9] z[6] z[9] + 
   256 Sqrt[2] a[7, 9] z[7] z[9] + 289 a[8, 9] z[8] z[9] + 
   196 Sqrt[6] a[4, 10] z[4] z[10] + 512 a[6, 10] z[6] z[10] + 
   289 Sqrt[3] a[7, 10] z[7] z[10] + 
   324 Sqrt[2] a[8, 10] z[8] z[10] + 361 a[9, 10] z[9] z[10] + 
   225 Sqrt[7] a[4, 11] z[4] z[11] + 
   289 Sqrt[5] a[6, 11] z[6] z[11] + 648 a[7, 11] z[7] z[11] + 
   361 Sqrt[3] a[8, 11] z[8] z[11] + 
   400 Sqrt[2] a[9, 11] z[9] z[11] + 441 a[10, 11] z[10] z[11] + 
   512 Sqrt[2] a[4, 12] z[4] z[12] + 
   324 Sqrt[6] a[6, 12] z[6] z[12] + 
   361 Sqrt[5] a[7, 12] z[7] z[12] + 800 a[8, 12] z[8] z[12] + 
   441 Sqrt[3] a[9, 12] z[9] z[12] + 
   484 Sqrt[2] a[10, 12] z[10] z[12] + 529 a[11, 12] z[11] z[12] +
    z[3] (81 Sqrt[3] a[3, 6] z[6] + 200 a[3, 7] z[7] + 
      121 Sqrt[5] a[3, 8] z[8] + 144 Sqrt[6] a[3, 9] z[9] + 
      169 Sqrt[7] a[3, 10] z[10] + 392 Sqrt[2] a[3, 11] z[11] + 
      675 a[3, 12] z[12]) + 
   z[5] (81 a[4, 5] z[4] + 121 a[5, 6] z[6] + 
      144 Sqrt[2] a[5, 7] z[7] + 169 Sqrt[3] a[5, 8] z[8] + 
      392 a[5, 9] z[9] + 225 Sqrt[5] a[5, 10] z[10] + 
      256 Sqrt[6] a[5, 11] z[11] + 289 Sqrt[7] a[5, 12] z[12])));

  vector z={z[1], z[2], z[3], z[4], z[5], z[6], z[7], z[8], z[9], z[10], z[11], z[12]};

The answer is very clear. I just want to know how to get the matrix a, the length-12 vector b, and scalar c in an efficient way.

|improve this question
$\endgroup$
  • $\begingroup$ Could you please post the relevant code you use to get the actual values for your {a,b,c}? $\endgroup$ – Dr. belisarius Mar 3 '13 at 8:30
6
$\begingroup$

One thing that may help is to make it a univariate function in a new variable t and extract derivatives. Below I changed your code to remove a space between vector and z on the last line or so.

With[{fexpt = fexp /. Thread[vectorz -> t*vectorz]},
 c = fexpt /. t -> 0;
 bvec = With[{der = D[fexpt, t] /. t -> 0}, Map[D[der, #] &, vectorz]];
 amat = With[{der = D[fexpt, {t, 2}]/2 /. t -> 0}, 
   1/2*Outer[D[der, ##] &, vectorz, vectorz]];
 {amat, bvec, c}]

I do not know if this is fastest possible, but it seems reasonable.

|improve this answer
$\endgroup$
  • $\begingroup$ This is indeed very fast, but it does not return correct values. $\endgroup$ – whuber Mar 3 '13 at 18:22
  • 1
    $\begingroup$ @whuber Actually it's a dearth of sharps, which by habit I have used only frugally since childhood. Should have been Outer[D[der, ##] &, vectorz, vectorz]]. Will edit. $\endgroup$ – Daniel Lichtblau Mar 5 '13 at 20:53
  • 1
    $\begingroup$ Strange. You're getting correct answers and your timing on this machine is 0.00388 seconds compared to 0.00943 on mine (based on 100 iterations of each). I'm still using MMA 8.0 on a four-core Xeon (about three years old). $\endgroup$ – whuber Mar 5 '13 at 21:37
  • 1
    $\begingroup$ @whuber I ran also using 100 iterations, on a Linux platform (that might or might not be relevant). In version 9 your's seems to be almost 3x faster. In version 8.0.4 yours remains unchanged but now they are nearly a dead heat. Maybe something got slower so I'll take a look. $\endgroup$ – Daniel Lichtblau Mar 5 '13 at 22:51
  • 1
    $\begingroup$ Okay, the slowdown is real. Reported as a bug. $\endgroup$ – Daniel Lichtblau Mar 5 '13 at 23:58
5
$\begingroup$

CoefficientRules seems to do a good job. Its results have to be transformed into actual arrays afterwards, but the timing is fast.

I am confident that those with greater facility in pattern manipulation can simplify this code, but it does work, albeit inelegantly. Most of the effort lies in constructing the matrix. The two forms of ij transmute expressions of the form {0,0,..,0,1,0,..,0,1,0,..0}->k to {{i,j}->k/2,{j,i}->k/2} and {0,...,0,2,0,...,0}->k to {i,i}->k, respectively, where i and j are the positions of the nonzero coefficients. Note the need in the first case to symmetrize the matrix by halving its off-diagonal values and replicating them across the diagonal. The output is cast into the form expected by SparseArray.

vector = Array[z, 12];
ij[Rule[x_, k_]] := With[{i = Flatten[Position[x, 1, 1, 2]]}, {i -> k/2, Reverse[i] -> k/2}];
ij[Rule[x : {___, 2, ___}, k_]] := With[{i = Position[x, 2, 1, 1]}, Flatten[{i, i}] -> k];
{a2, b1, c0} = Block[{r, a2, b1, c0},
    r = CoefficientRules[fexp, vector];
    c0 = First[Select[r, Plus @@ (#[[1]]) == 0 &] /. Rule[_, c_] :> c];
    b1 = Plus @@ (Select[r, Plus @@ (#[[1]]) == 1 &] /. 
        Rule[i_, x_] :> x i);
    a2 = Select[r, Plus @@ (#[[1]]) == 2 &];
    a2 = SparseArray[Flatten[a2 /. x_Rule :> ij[x]]];
    {a2, b1, c0}
    ]; // AbsoluteTiming

$\{0.0090005, Null\}$

To check, let's compare the value calculated with {a2, b1, c0} to fexp:

vector . a2 . vector + b1.vector + c0 - fexp // Simplify

$0$

It truly is the same.

This method assumes the input really is a quadratic form in the z[..] (including linear forms and constants as special cases). If not, the output will be incorrect.

|improve this answer
$\endgroup$
  • $\begingroup$ Thanks! I learn a lot from your code. $\endgroup$ – Orders Mar 6 '13 at 2:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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