How to create matched "minimum phase" for a system of parallel FIRs? Problem statement

I have a collection of magnitude (only) responses I'd like to turn into FIR filter kernels that are


matched in phase
"minimal" in phase, with respect to the complete system
have little to no "pre-ring"
maintain the design magnitude to a high degree


In the ideal case, I'm hoping for a solution that is analytic rather than parametric.

Meeting the first requirement can be expressed as:

$$
\phi_s(f) = arg\{ H_0(f) \} = arg\{ H_1(f) \} = ... = arg\{ H_n(f) \}
$$

And, is very easy to meet in a linear phase case.

Accumulating phase

In attempting to meet the second condition, I've tried something akin to a "cascading" of phase given minimum phase versions of each filter. I've done this by


returning the "Hilbert" minimum phase for each magnitude response
accumulating these responses, and then defining this as the system phase
reassigning system phase to each response


Here's an expression for the accumulation I used to find the system phase:

$$
\phi_s(f) = \sum_{i=0}^n - \mathscr H \{ \ln( |H_n(f)| ) \}
$$

Unfortunately, for my responses, this results in a poor outcome, returning kernels with way too much smearing / lag at low frequencies in the time domain. All this smearing is far from "minimal" in any sense!

I'm suspicious this may be a result of the two definitions of minimum phase discussed earlier.

Minimum phase of L-2 norm of magnitudes

The most promising results I've had after a bit of hacking is to take the L-2 norm of my magnitude set and use this as a minimum phase target:

$$
|H_{L2}(f)| = \biggl(\sum_{i=0}^n |H_n(f)|^2 \biggr)^{1/2}
$$

$$
\phi_{L2}(f) = - \mathscr H \{ \ln (|H_{L2}(f)|) \}
$$

Setting $\phi_s(f) = \phi_{L2}(f)$ returns kernels with coefficients packed towards the beginning and end of the table. Rotating the table, a linear phase offset, by $\max\bigl(\phi_{L2}(f)\bigr)$, turns the end into a "pre-ring".

While not beautiful, this does get me closer to what I think I'm looking for.

Some more details on my magnitudes

$|H_0(f)|$ is allpass; the remaining are highpass with every increasing roll-off slopes and cut-off frequencies.

Questions

Is there a nice off the shelf know solution? In my mind, the problem is very similar to a matched phase cross-over.

Am I asking too much? In the ideal case I'd like to have no "pre-ring". Maybe that's not possible?

Am I asking the right question(s)? In other words, is "minimum phase" what I'm looking for. What I really want is for my kernels to be matched in phase, and respond "as soon as possible".
The Hilbert "transform" or relation define one phase response for a given magnitude response, so you can't get both matched and minimum phase in your case.

Regards.