Context: I am creating a special ballistics simulation with nonstandard physics, and I have reached a step where I need to find the global minimum of a rational function very fast to be able to run the simulation interactively in real-time. I have documented my research and implementation in C#: it can be found on GitHub.
Goal: Find the global minimum of a univariate, nonnegative rational function of the form $$\left\| v(T)\right\|^2:= \left\| \frac{x(T)}{T} \right\|^2=\frac{x(T)\cdot x(T)}{T^2}$$
for some Taylor polynomial with vector "coefficients" $$x(T):=\sum_{k=0}^{n} \frac{T^k}{k!}\vec{a_k}$$ with $x(T)$ and therefore $a_k$ of arbitrary dimension.
I am interested in the domain of $\mathbb{R}_{>0}$. As you notice, this function is strictly nonnegative on this interval.
My first thought, which is a standard approach for optimization, is to use calculus to take the derivative and set it to $0$ to find critical points and test for the minimum (the function is differentiable on $\mathbb{R}_{>0}$ so the critical points must be where the derivative is $0$). So I obtained: $$2\frac{x(T)}{T}\cdot \frac{x(T)-T \frac{dx}{dT}(T)}{T^2} = 0$$ Based on that, I defined a new function $\text{Critical}$ whose roots constitute critical points of the original function: $$\text{Critical}(T):=x(T)\cdot \left(T \frac{dx}{dT}(T)-x(T)\right) = 0$$ This is a polynomial in $T$. However, I quickly encountered a problem: how do I find all the positive real roots of the derivative polynomial reliably? I searched for algorithms and found some.
However, I then rapidly hit another obstacle: many of these root-finding algorithms rely on operations on the polynomial's coefficients. This meant that I was forced to expand the derivative to find its coefficients. Obtaining that the coefficient of $T^k$ is $$\sum_{j=0}^k \frac{\vec a_j \cdot \vec a_{k-j}(1-k+j)}{j!(k-j)!}$$
where $\vec a_j \cdot \vec a_{k-j}$ is a dot product for $k=0,1,...,n$, or in other words, for all the vectors in defining the Taylor polynomial of $x(T)$.
As you can notice, the time complexity of this algorithm is $O(n^2)$ dot products, which is quite expensive computationally. On top of that, I need to apply a root-finding algorithm on this new polynomial to find all its roots then test each critical point in the original function, which is not cheap either, resulting in a lot of latency. The current root-finding algorithm that I am using is the continued fractions algorithm based on Vincent's theorem and Descartes' rule of signs. It's an iterative interval bracketing method relying on bounds on the positive real roots of polynomial to initialize the interval. More on it here and you can find my C# implementation here. I've tried other iterative methods such as Newton's but they were not stable enough to guarantee convergence.
Also, the accumulation of errors from the evaluation of the polynomial adds up to quite a lot, making it prone to numerical instability. It has to pass through the many computations of the Taylor polynomial evaluation with the exponentiation and vector scalings, then the dot product, etc. so it would be preferable to reduce the amount of steps and/or evaluations needed.
Upon investigation, I have found a technique called Semi-definite optimization (also relevant: here) that could be used to minimize rational functions directly through representation of a nonnegative polynomial as a sum of squares. Nevertheless, this also seems quite involved, having to compute a matrix describing an affine space on a monomial basis, then making sure that each term is nonnegative. I am not familiar with this, and if anyone could suggest if this is potentially a better solution to this problem, then please feel free to explain.
After some further analysis, if I understand this properly, it seems like you could find the global minimum of a rational function by reformulating it into a Semi-definite programming problem, but the specifics are not mentioned, and I do not know how you do that. If someone is more familiar with this topic and could give advice on its effectiveness, it would be much appreciated. I believe directly aiming to find the global minimum should be faster to avoid having to find all the roots of the derivative of the function.
I am looking for suggestions on optimization in terms of performance. If anyone could offer ideas to be able to minimize the original function fast, it would be greatly appreciated.
As a side note, following this problem, I would like generalize this problem to any function of the form $$x^{(k)}(T):=\frac{x(T)\cdot x(T)}{T^{2k}}$$
I am guessing that if the original problem is solved efficiently, it should not be too difficult to extend the logic to this one, but this comes after solving the first problem in an efficient manner.
There are many methods for numerical optimization. I would suggest you experiment with them. One standard approach is gradient descent. Another is Newton's method.
Let $$f(T) = \left\lVert{x(T) \over T}\right\rVert^2.$$ Gradient descent requires the ability to evaluate $f(T)$ and to evaluate the first derivative $f'(T)$. Notice that we can obtain an algebraic expression for $f'(T)$: $$f(T) = \sum_j \left({x_j(T) \over T}\right)^2,$$ so $$\begin{align*}f'(T) &= \sum_j 2 {x_j(T) \over T} \left({x'_j(T) \over T} - {x_j(T) \over T^2}\right)\\ &= 2 {x(T) \cdot x'(T) \over T^2} - 2 {\lVert x(T)\rVert^2 \over T^3}, \end{align*}$$
where $$x'(T) = \sum_{k=1}^n {T^{k-1} \over (k-1)!} \overline{a}_k.$$
Therefore, you can evaluate $f(T)$ and $f'(T)$ using $2n$ vector-sums and $3$ dot-products. This is enough to use gradient descent.
If you expand the calculation a bit further, you can compute the second derivative $f''(x)$ and then use Newton's method. The advantage of Newton's method is that it converges much faster, often in only a few iterations. I suspect this will make up for the additional expense of computing the second derivative.