Shape optimization has become a key tool in structural design, enabling the development of efficient and lightweight structures, particularly in engineering applications. The shape derivative plays a crucial role in shape optimization by providing sensitivity information of an objective function with respect to infinitesimal perturbations of the structural boundary. In recent years, shape derivatives have been successfully applied to a wide range of engineering problems, often in combination with level-set methods driven by Hamilton–Jacobi advection or geometric optimization techniques. When dealing with constrained optimization problems, however, the shape derivative must be embedded within a global numerical scheme capable of handling many constraints. A commonly adopted approach is the augmented Lagrangian method, which may suffer from slow convergence and requires careful parameter tuning. More recently, shape derivatives have been incorporated into the null-space algorithm, which allow for the efficient treatment of multiple constraints with a reduced number of parameters. In its standard formulation, null-space relies on a linearization of both the objective functional and the constraints. Such approach has shown promising performance when combined with level-set approaches and Hamilton–Jacobi-based updates. Motivated by these developments, we propose an extension of the null-space algorithm for optimization problems in which functionals involving the solution of partial differential equations are linearized, while purely geometric functionals are treated in a nonlinear fashion, taking advantage of their lower computational cost. The potential efficiency of the proposed approach is explored through a series of numerical examples.