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.