TaoNullSpace.hpp

Go to the documentation of this file.

/** \file TaoNullSpace.hpp
 * \brief TAO shell interface for null-space constrained optimization.
 */
 
#ifndef MOFEM_TAO_NULL_SPACE_HPP
#define MOFEM_TAO_NULL_SPACE_HPP
 
namespace MoFEM {
 
/** Registered TAOSHELL-based null-space constrained optimizer.
 * \ingroup mofem_petsc_solvers
 */
inline constexpr char TAONULLSPACE[] = "nullspace";
 
/** Register TAONULLSPACE before TaoSetType() or TaoSetFromOptions().
 *
 * The shell consumes the same objective/gradient, inequality and inequality
 * Jacobian callbacks as TAOALMM, with c_i >= 0. Both combined and separate
 * objective/gradient callbacks are supported. The null-space flow follows
 * Feppon, Allaire and Dapogny, ESAIM: COCV 26 (2020), article 90.
 * Objective-only trial evaluations use the parent objective and constraint
 * values without requesting the objective gradient or constraint Jacobian.
 * Applications providing only a combined callback retain PETSc's combined
 * evaluation fallback. The default backtracking and the More--Thuente domain
 * safeguard request a derivative only after a trial passes the merit check;
 * native More--Thuente still evaluates the derivative for its curvature test.
 *
 * At each accepted iterate, rows with c_i <= h * ||J_i||_2 form the multiplier
 * candidates, following Section 4.1, equation (4.3), of Feppon et al. Here
 * h = active_tol is a control-space activation distance in the Euclidean Vec
 * metric. This replaces the former uniform tolerance on constraint values.
 * Positive scaling of a constraint and its Jacobian row leaves this candidate
 * classification unchanged. A nonnegative least-squares solve fits
 * J_c^T lambda to grad f; multipliers above multiplier_tol select the binding
 * rows B. The restoration rows R contain B and all rows with c_i <= 0.
 * Ordered row selections and their accepted-iterate Jacobians are copied
 * before constructing C C^T inverse
 * actions for projection and restoration.
 *
 * The default inner TAOSHELL takes one flow step in direction -grad Phi at
 * the accepted base point, with projection and merit data frozen. The default
 * "nullspace" line search follows Algorithm 1: halve the step until the
 * frozen merit strictly decreases, or accept the final valid trial when the
 * trial budget is reached. That terminal fallback may increase the merit.
 * Invalid physical evaluations cannot be accepted. After the inner solve
 * returns, the accepted outer update immediately rebuilds row selections and
 * multipliers, including constraints newly violated by the step. Convergence uses
 * stationarity, complementarity and primal feasibility through standard TAO
 * monitoring and convergence callbacks.
 * An iterative multiplier solve reaching its iteration limit emits a warning
 * and uses its current estimate; the outer convergence checks still apply. Other
 * multiplier-solver failures remain errors.
 *
 * Options following the parent TAO prefix:
 * - -tao_nullspace_monitor_kkt: print stationarity, projected-gradient and
 *   complementarity Euclidean norms at each outer iteration, default false.
 * - -tao_nullspace_alpha_j: objective scale, default 1.
 * - -tao_nullspace_alpha_c: restoration scale, default 1.
 * - -tao_nullspace_active_tol: finite nonnegative control-space activation
 *   distance h, default 1e-6; h = 0 selects saturated or violated candidates.
 * - -tao_nullspace_multiplier_tol: binding-row tolerance, default 1e-10.
 * - -tao_nullspace_min_step: minimum accepted step norm, default 1e-14.
 *
 * The inner solver uses the additional prefix "tao_nullspace_". Its default
 * TAO type is "shell" and its line-search type is "nullspace".
 * -tao_nullspace_tao_ls_stepinit sets the initial time step (default 1), and
 * -tao_nullspace_tao_ls_max_funcs sets the trial budget (PETSc default 30).
 * The last valid trial may return TAOLINESEARCH_SUCCESS_USER without a merit
 * decrease; the outer solver honors this fallback. Invalid/nonfinite physical
 * data and a step below the configured minimum still cause failure.
 * A completed shell step returns TAO_CONVERGED_USER; the outer TAO checks KKT.
 * Explicit LMVM or other line-search choices are alternatives to this default
 * Algorithm 1 step. Iterative inner types honor -tao_nullspace_tao_max_it.
 * LMVM history recycling follows -tao_nullspace_tao_lmvm_recycle (PETSc
 * default false). With one inner iteration, recycling allows curvature history
 * to accumulate across outer iterations. The frozen merit changes when its
 * active sets, multipliers or restoration matrix change, so retained history
 * is an approximation; the outer acceptance and KKT checks still apply.
 * With More--Thuente, invalid physical data stop cubic interpolation and
 * trigger geometric backtracking using the true merit and gradient.
 * Its sufficient-decrease safeguard does not guarantee Wolfe curvature.
 * Every option above also receives the parent TAO prefix.
 *
 * Multiplier selection defaults to TAOBLMVM under "tao_nullspace_dual_".
 * For orthogonal nonzero constraint rows, select
 * -tao_nullspace_dual_tao_type diagonal. TAODIAGONALNNLS then evaluates
 * lambda = max(0, (J * grad f) / diag(J J^T)) directly, without an iterative
 * solve or line search, as in Remark 3.10 of Feppon et al. Cellwise P0 plastic
 * rows have disjoint control supports, which satisfies this condition even
 * after volume scaling. The diagonal solver checks the normalized Gram matrix
 * for orthogonality; coupled rows require a general solver such as TAOBLMVM.
 * Its monitor and convergence test use the actual bound-projected gradient.
 * The reusable API is declared in TaoDiagonalNNLS.hpp:
 * TaoRegisterDiagonalNNLS() and TaoDiagonalNNLSSetProblem(tao, J, target)
 * define min_{lambda >= 0} 0.5 * ||J^T lambda - target||^2 for type
 * TAODIAGONALNNLS ("diagonal").
 * For a P0 constraint c = sqrt(V) * (delta_kappa - alpha), the Jacobian-row
 * norm is sqrt(V) * sqrt(1 + ||grad alpha||_2^2), so volume cancels from candidate
 * selection. This does not normalize the multiplier, feasibility or stopping
 * tolerances.
 *
 * Projection KSPs use
 * "tao_nullspace_binding_" and "tao_nullspace_restoration_"; both default to
 * PREONLY with LU factorization through MUMPS. These prefixes all follow the
 * parent prefix, and PETSc options can override the defaults.
 * The MPI regression sets mat_mumps_icntl_20 to 0 under both KSP prefixes
 * for centralized RHS input with empty local active-row partitions.
 *
 * This implementation handles the supplied inequality system; separate TAO
 * variable bounds and equality constraints are rejected. Express
 * variable bounds as inequality rows. The nonlinear epigraph regression
 * exercises smooth points. At a cone apex, choosing an appropriate generalized
 * Jacobian remains the application's responsibility and is not covered by that
 * test. Selected-row Gram systems must admit the configured inverse action.
 */
MoFEMErrorCode TaoRegisterNullSpace();
 
/** Bind the inequality storage already registered on the parent TAO.
 *
 * Call after selecting TAONULLSPACE and before TaoSetUp(). The shell retains
 * PETSc references to all supplied objects, and releases them on a type change
 * or parent destruction. The constraint DM carries application numbering and
 * may be null for purely algebraic problems. A null preconditioner uses the
 * Jacobian. Callbacks remain registered through the standard TAO setters.
 */
MoFEMErrorCode TaoNullSpaceSetConstraints(Tao tao, DM constraint_dm,
                                          Vec constraints, Mat jacobian,
                                          Mat preconditioner = nullptr);
 
/** Borrow the owned inner TAO for configuration; do not destroy it.
 * After parent TaoSetUp(), it has its own solution and gradient vectors and
 * separate and combined objective/gradient callbacks for the frozen merit
 *
 * Phi(x) = alpha_j [f(x) - lambda_B^T c_B(x)]
 *          + alpha_c/2 c_R(x)^T S_R c_R(x).
 *
 * Here S_R = (C_R C_R^T)^(-1), and lambda_B is the projection multiplier
 * obtained from the accepted-iterate binding Jacobian and objective gradient.
 * These quantities remain fixed while trial f and c are evaluated. Trial
 * gradients and Jacobians are evaluated only for gradient or combined
 * callbacks; the inner gradient is the derivative of Phi, not a separately
 * projected trial gradient.
 *
 * Callbacks require a frozen outer iterate established by TaoSolve(); calling
 * them earlier returns PETSC_ERR_ARG_WRONGSTATE. Binding/restoration rows and
 * projection data remain frozen until the inner solve returns. The default
 * shell takes exactly one flow step. An explicitly selected iterative solver
 * resets history between outer iterations and may use it within its configured
 * inner solve. A retained inner TAO loses
 * these callbacks when the parent context is destroyed or changes type.
 */
MoFEMErrorCode TaoNullSpaceGetSubsolver(Tao tao, Tao *subsolver);
 
/** Borrow the inequality multiplier vector after TaoSetUp().
 * It is initially zero and uses the parent constraint layout and the
 * nonnegative-multiplier convention for c_i >= 0. Applications may initialise
 * it as for TAOALMM; the outer solver updates it through nonnegative
 * least-squares selection at each accepted iterate. The stationarity balance
 * is grad f - J_c^T lambda = 0, so scaling a constraint also changes the units
 * of its returned multiplier.
 */
MoFEMErrorCode TaoNullSpaceGetMultipliers(Tao tao, Vec *multipliers);
 
} // namespace MoFEM
 
#endif // MOFEM_TAO_NULL_SPACE_HPP