Factorising blinear forms with Mathematica

Question

Say I have some symbolic output that looks like the following:

expr = A a1 b1 + B a1 b2 + C a2 b1 + D a2 b2

I would like to factorise this in the following way:

Input: FactorQuadratic[expr, {a1, a2}, {b1, b2}]

Output: { { A, B }, { C, D } }

In such a way that $\begin{bmatrix} a_1 & a_2\end{bmatrix} \begin{bmatrix} A & B \\ C & D \end{bmatrix}\begin{bmatrix}b1\\b2\end{bmatrix}$ is equal to the original expression.

Is there a function that does this, or something like it? The closest I'm aware of is FactorList, but I'm not entirely sure how to slice and dice indices on it to make it do what I want.

Answer 1

Per this answer, you can compute the Hessian:

D[A a1 b1 + B a1 b2 + C a2 b1 + D a2 b2, {{a1, a2}}, {{b1, b2}}]
(* {{A, B}, {C, D}} *)

Answer 2

Maybe this:

expr = $A a1 b1 + $B a1 b2 + $C a2 b1 + $D a2 b2;
v1 = {a1, a2};
v2 = {b1, b2};

Last@CoefficientArrays[#, v2] & /@ 
  Last@CoefficientArrays[expr, v1] // Normal

(*  {{$A, $B}, {$C, $D}}  *)

Note C and D are protected symbols and not meant to be used as user parameters.