<think> We are given a complex optimization problem involving a mesh with constraints to prevent triangle inversion. The user has a function to compute the maximum step size that ensures validity (no inverted triangles) and wants to use this in Mathematica's `FindMinimum` with conjugate gradient method.

The key issue is that the built-in line search in FindMinimum does not seem to accept a user-defined maximum step size in the way the user attempted. The user tried using a delayed rule (:>) for "MaxRelativeStepSize" but it resulted in an error.

After reviewing Mathematica's documentation, the official methods for FindMinimum do not expose a way to inject a custom step size control that depends on the current point. However, there might be alternative approaches:

1. Implementing a custom line search procedure: This would involve more work but gives full control. The user would have to implement the conjugate gradient (or other method) manually, using their own line search that respects the maximum step size.

2. Using NDSolve for gradient flow: Another idea is to set up a differential equation that follows the gradient, and use WhenEvent or similar to enforce constraints. However, this might be less efficient and not directly applicable.

3. Using the "TrustRegion" method: The trust region method inherently restricts steps to a region. If we can define the region by the maximum step size, that might work. However, the trust region is typically a sphere and the user's constraint is more complex (depending on the geometry of the mesh). But note: the user's constraint is that for a given current point X, the step must be less than S (a scalar) to remain valid. This is a norm constraint? Actually, the user states: "calculate a maximum steps-size S such that for any actual step size s (0 < s < S) the resulting vertex coordinates will also be valid". This suggests that the step in the direction of the gradient (or the search direction) must be bounded by S to avoid inversion. This is a one-dimensional constraint along the line.

Therefore, the line search must be constrained to step sizes in the interval (0, S). The built-in line search in Mathematica does not allow for such a constraint.

4. Using FindMinimum with StepMonitor or EvaluationMonitor: We might try to interrupt the process and adjust the step size, but this is not straightforward and might break the algorithm.

Given the above, the most promising approach is to implement a custom gradient descent (or conjugate gradient) method. Since the user already has the gradient and the potential function, and a way to compute the maximum step size S for the current point and direction, they can build their own optimizer.

However, the user specifically asked for a way to do this within FindMinimum. Let me check for any undocumented options.

After some research, I recall that the "LineSearch" method has an option "MaxRelativeStepSize" that can be set to a constant. But the user wants to compute it dynamically at each step. The error message suggests that the delayed rule is not being accepted.

Looking at the documentation for "LineSearch" in the FindMinimum context:

• "MaxRelativeStepSize" is an option for "LineSearch" that specifies the maximum size of the step relative to the current point. It is a machine number, so it must be a number, not a function.

Therefore, we cannot use a function to compute it at each step.

Alternative idea: use the "TrustRegion" method and set the trust region radius to the maximum step size? But note:

• The trust region method uses a region (usually Euclidean) and solves a subproblem. The user's constraint is not a Euclidean ball but a condition that the step must be less than S to avoid inversion. However, if we set the trust region radius to S (the maximum step that is safe) then the trust region step would be bounded by S. But note: the trust region step is computed by approximately minimizing a quadratic model within the trust region. This might not be the same as the line search, but it does bound the step.

But the user's constraint is that for the specific direction of the gradient (or conjugate direction), the step must be less than S. The trust region uses a spherical constraint, so it might take a step in a different direction? Actually, the trust region step is not necessarily along the gradient. So this might not be equivalent.

Moreover, the trust region radius is usually adjusted automatically, and we would need to set it at each step to the computed S. But the trust region method in FindMinimum does not allow the radius to be set by a function.

Given these limitations, it seems that the built-in methods do not support the required dynamic step size constraint.

Therefore, I propose a custom approach:

The user can implement a simple gradient descent with adaptive step size, using their CalculateMaxRelativeStepSize function. The algorithm would be:

1. Start at X0.
2. Compute the gradient g at current X.
3. Compute the search direction d (for gradient descent, d = -g; for conjugate gradient, we would need to update the direction).
4. Compute the maximum allowed step size S for the current X and direction d.
5. Perform a line search in the direction d for a step size s in (0, S] that minimizes the potential (or satisfies Wolfe conditions, etc.). Since the potential is expensive, we might use a backtracking or interpolation-based line search.

The user already has a function to compute S, so the main challenge is the line search. However, the user might not need a full line search: they could use a fixed fraction of S and then adjust if necessary? But that might be inefficient.

Alternatively, the user could use a simple backtracking line search that starts with s = min(S, s0) (where s0 might be 1 or a value from previous steps) and then reduces until the potential decreases (and the step remains valid). But note: the user's S already ensures validity, so we only need to ensure the Armijo condition (sufficient decrease) for the step.

The algorithm for one iteration:

• Let α be a reduction factor (e.g., 0.5) and c a small constant (e.g., 10^-4) for the Armijo condition.
• Set s = min(S, initial_step) (initial_step could be 1 or based on previous steps).
• Evaluate f_new = potentialFn(X + s * d).
• While f_new > f_current + c * s * g.d (i.e., not sufficient decrease) and s is not too small:
  • Reduce s by s = α * s
  • Re-evaluate f_new.
• Then update X = X + s * d.

But note: the user's S ensures that the step s in the direction d will not invert triangles. So we are safe as long as s <= S.

For conjugate gradient, we would need to store the previous gradient and update the direction. The Polak-Ribiere or Fletcher-Reeves formulas can be used.

Given the complexity of the problem and the need for efficiency (25,000 vertices), a well-implemented conjugate gradient with custom line search might be necessary.

But the user asked for a way within FindMinimum. Since that seems impossible, I must conclude that a custom implementation is the way to go.

However, let me check if there is a way to use FindMinimum with a user-defined line search. The "StepControl" option only allows built-in methods. There is no hook for a user-defined line search.

Another idea: use FindMinimum with Method -> "Gradient" and control the step size via StepMonitor? But we cannot override the step size from StepMonitor.

Alternatively, we can use FindMinimum in a step-by-step manner by using "StepMonitor" to record the steps and then manually adjust the next point? But that would break the algorithm.

Given the above, I think the best solution is to implement a custom optimizer.

But to help the user, I can sketch a simple gradient descent with adaptive step:

(* Pseudocode *)

customGradientDescent[potentialFn_, gradientFn_, maxStepFn_, X0_, maxIter_: 100, tol_: 10^-6] := 
Module[{X, g, f, fOld, s, S, d, iter = 0, α = 0.5, c = 0.0001, initialStep = 1.0},
  X = X0;
  f = potentialFn[X];
  g = gradientFn[X];
  While[iter < maxIter && Norm[g] > tol,
    iter++;
    d = -g; (* gradient descent direction *)
    S = maxStepFn[X, d]; (* maximum step to avoid inversion *)
    s = Min[initialStep, S];
    (* Backtracking line search *)
    While[True,
      Xnew = X + s * d;
      fNew = potentialFn[Xnew];
      If[fNew <= f + c * s * (g.d), 
        Break[]; (* sufficient decrease *)
      ];
      s = α * s;
      If[s < 10^-10, Break[];]; (* avoid infinite loop *)
    ];
    (* Update *)
    X = Xnew;
    fOld = f;
    f = fNew;
    g = gradientFn[X];
    (* Adjust initialStep for next iteration? Maybe set to 2*s or something *)
    initialStep = 2*s; (* or keep as 1? *)
  ];
  {X, f}
]

This is a very basic gradient descent. For conjugate gradient, we would need to store the previous gradient and update the direction. Here is a sketch for conjugate gradient:

customConjugateGradient[potentialFn_, gradientFn_, maxStepFn_, X0_, maxIter_: 100, tol_: 10^-6] := 
Module[{X, g, gOld, f, d, s, S, iter = 0, α = 0.5, c = 0.0001, initialStep = 1.0, β},
  X = X0;
  g = gradientFn[X];
  d = -g;
  f = potentialFn[X];
  While[iter < maxIter && Norm[g] > tol,
    iter++;
    S = maxStepFn[X, d];
    s = Min[initialStep, S];
    (* Backtracking line search *)
    While[True,
      Xnew = X + s * d;
      fNew = potentialFn[Xnew];
      If[fNew <= f + c * s * (g.d), 
        Break[];
      ];
      s = α * s;
      If[s < 10^-10, Break[];];
    ];
    X = Xnew;
    gOld = g;
    g = gradientFn[X];
    (* Fletcher-Reeves: *)
    β = (g.g) / (gOld.gOld);
    (* Or Polak-Ribiere: β = Max[0, (g.(g - gOld)) / (gOld.gOld)]; *)
    d = -g + β * d;
    (* Check if we need to restart: for example, every n steps or if not descent direction *)
    If[d.g >= 0, (* not a descent direction, restart *)
      d = -g;
    ];
    initialStep = 2*s; (* or adjust based on previous step *)
  ];
  {X, f}
]

Note: The conjugate gradient method should restart periodically (e.g., every n steps, n=dimension) but in high dimension we might restart every n iterations. Also, the step size adjustment for initialStep is heuristic.

The user must provide:

• potentialFn[X]: computes potential at X.
• gradientFn[X]: computes gradient at X.
• maxStepFn[X, d]: computes the maximum step size S for the current X and direction d such that for any s in (0,S] the new point X+s*d is valid.

This custom approach should work and allow the user to enforce the constraint.

In summary, since FindMinimum does not support a dynamically computed maximum step size in the line search, the user must implement their own optimization routine. The above provides a template for gradient descent and conjugate gradient methods.

I hope this helps. </think>

Comprehensive Answer

The core issue is that Mathematica's built-in FindMinimum function does not support dynamically computed maximum step sizes during line searches, which are essential for your mesh optimization to prevent triangle inversion. Your attempt to use "MaxRelativeStepSize" :> CalculateMaxRelativeStepSize[...] fails because "MaxRelativeStepSize" requires a static numeric value, not a function that updates at each step. Below, I outline a structured solution and alternatives.

---

1. Why the Built-in Approach Fails

• "MaxRelativeStepSize" Limitation: This option only accepts a fixed numeric value (e.g., 0.5). It cannot evaluate a function like CalculateMaxRelativeStepSize at each iteration.
• Error Analysis: The error FindMinimum::lsmrss occurs because FindMinimum expects a number, not a delayed rule (:>) or function.
• No Direct Hook in Line Search: Mathematica's line-search methods (e.g., "LineSearch" in "ConjugateGradient") lack callbacks for user-defined step-size constraints.

---

2. Workaround: Custom Optimization Algorithm

Since FindMinimum cannot enforce dynamic step-size constraints, implement a custom gradient-based optimizer (e.g., Gradient Descent or Conjugate Gradient) with your CalculateMaxRelativeStepSize function. This approach guarantees validity and leverages your existing gradient computations.

Pseudocode for Custom Conjugate Gradient

CustomConjugateGradient[potentialFn_, gradientFn_, maxStepFn_, X0_, maxIter_: 100, tol_: 1*^-6] := 
Module[{
  X = X0, 
  g = gradientFn[X0], 
  d = -gradientFn[X0], (* Initial direction: negative gradient *)
  f = potentialFn[X0],
  iter = 0, 
  α = 0.5, (* Backtracking reduction factor *)
  c = 10^-4, (* Armijo condition constant *)
  s, S, fNew, XNew, gOld, β
},
  While[iter < maxIter && Norm[g] > tol,
    iter++;
    S = maxStepFn[X, d]; (* Max step to avoid inversion *)
    s = Min[1.0, S]; (* Start with s <= S *)
    
    (* Backtracking line search: enforce decrease & validity *)
    While[True,
      XNew = X + s*d;
      fNew = potentialFn[XNew];
      If[fNew <= f + c*s*(g.d), 
        Break[]; (* Armijo condition met *)
      ];
      s *= α; (* Reduce step size *)
      If[s < 10^-10, Break[]]; (* Avoid infinite loop *)
    ];
    
    (* Update position and gradient *)
    X = XNew;
    gOld = g;
    g = gradientFn[X];
    
    (* Fletcher-Reeves conjugate direction update *)
    β = (g.g)/(gOld.gOld); 
    d = -g + β*d;
    
    (* Optional: Reset to steepest descent if not descent direction *)
    If[d.g >= 0, d = -g];
  ];
  {X, f} (* Return solution *)
]

Key Components:

1. maxStepFn: Your existing CalculateMaxRelativeStepSize function, which ensures no triangle inversion for steps s < S.
2. Backtracking Line Search: Starts with s = Min[1.0, S] and reduces s until the Armijo condition (sufficient function decrease) is met.
3. Conjugate Direction Update: Uses the Fletcher-Reeves formula. Restarts with steepest descent if d is not a descent direction.
4. Validity Guarantee: Since s ≤ S, all intermediate points are valid (no inverted triangles).

---

3. Integration with Your Problem

• Functions to Define:
  • potentialFn[X_]: Computes potential from vertex coordinates X.
  • gradientFn[X_]: Computes analytical gradient.
  • maxStepFn[X_, d_]: Computes max step S for direction d at X (prevents inversion).
• Initialization:

  {XOpt, fOpt} = CustomConjugateGradient[potentialFn, gradientFn, maxStepFn, X0, 100];
  

• Efficiency: Warm-start s using values from prior iterations (e.g., s = 2*previous_s).

---

4. Why Not Use Built-in Methods?

• "TrustRegion": Restricts steps to a Euclidean ball, but your constraint is directional (along d). Trust regions may overshoot S or compute inefficient steps.
• "Gradient" Method: Allows fixed "StepSize" but not dynamic adjustments. Impractical for your adaptive needs.
• StepMonitor/EvaluationMonitor: Can log data but cannot modify step sizes mid-iteration.

---

5. Validation and Tuning

• Test on Small Mesh: Verify correctness with a toy mesh (e.g., 10 vertices).
• Convergence Criteria: Use Norm[g] < tol or relative change in X.
• Parameter Tuning:
  • α (backtracking factor): Typical values 0.1–0.8.
  • c (Armijo constant): Typical values 10^-4–10^-2.
  • maxIter: Set based on mesh size (e.g., 100–1000).

---

6. Alternatives (If Custom Code Is Infeasible)

• Fixed Step Size: Use FindMinimum with Method -> {"Gradient", "StepSize" -> s_fixed}, where s_fixed is a conservative fraction of the minimum S observed in test runs. Drawback: Slow convergence.
• External Libraries: Link to high-performance optimizers (e.g., IPOPT via MATLink) that support custom constraints. Overkill for most cases.

---

Conclusion

Implement a custom Conjugate Gradient optimizer with backtracking line search and your CalculateMaxRelativeStepSize function. This is the only robust way to enforce dynamic step-size constraints and prevent triangle inversion in Mathematica. Example pseudocode is provided above for direct adaptation.

For large meshes (~25k vertices), optimize performance by:

• Using compiled functions for potentialFn, gradientFn, and maxStepFn.
• Parallelizing gradient/potential computations over mesh partitions.
• Storing edge lengths/angles incrementally to avoid redundant calculations.