Practically, minimizing the objective function \(Q_{\lambda }(\psi )\)  with respect to \(\psi \)  poses several interconnected challenges: the loss function \(\ell \)  is non-convex due to the marginalized random frailty; the penalty functions \(p_{\lambda }\)  may also be non-convex; and, the fusion penalty component does not admit standard algorithms for general fused Lasso tailored to linear regression [1]}. Finally, the combination of penalties requires tuning a two-dimensional regularization parameter \(\lambda \) . In this section, we propose a comprehensive optimization routine to simultaneously and efficiently handle these challenges.
