Polynomial optimization, how to solve numerically?

Question

For context this problem is derivative from this problem regarding polynomial position estimation.

---

I have a problem where I have two vectors , $\bf v$ and $\bf w$.

Let us assume that they are linked like so:

$${\bf w} = {\bf Mv} = \left[\begin{array}{ccccc} 1 &&&&\\a&1\\a^2&2a&1\\a^3&3a^2&3a&1\\\vdots&\ddots&&\ddots&1\end{array}\right] {\bf v}$$ Now I want to find $a$. It is possible that no exact solution exists, so maybe I will have to estimate $a$ as$$\min_{a}\|{\bf Mv-w }\|$$

We can assume for simplicity that $a \in [-1,1]$.

How can this be done?

---

Own work As the search space is one dimensional real variable $a$ living in a confined space. The most obvious and simple approach I can think of would be to try a binary interval split. But I suspect more clever ways can be used which utilize the behavior of this problem.

Answer 1

Calling $\mathcal{O}(a) = \|w-M(a)\cdot v\|^2$ we can use a simple gradient steepest descent method do locate $\min_a \mathcal{O}(a)$. Choosing $\mu$ small we can follow with

$$ a_1 = a_0 -\mu \nabla_a \mathcal{O}(a_0) $$

etc.

Follows a python script showing the procedure for $\text{rank}\left(M(a)\right) = 5$

a0 = -0.8

import numpy as np
import matplotlib.pyplot as plt

def A(a):
    a2 = a*a
    a3 = a*a2
    a4 = a*a3
    return np.array([[1,    0,  0,   0,0],
                     [a,    1,  0,   0,0],
                     [a2, 2*a,  1,   0,0],
                     [a3,3*a2,3*a,   1,0],
                     [a4,4*a3,6*a2,4*a,1]])

def dA(a):
    a2 = a*a
    a3 = a*a2
    return np.array([[   0,    0,   0,0,0],
                     [   1,    0,   0,0,0],
                     [ 2*a,    2,   0,0,0],
                     [3*a2,  6*a,   3,0,0],
                     [4*a3,12*a2,12*a,4,0]])   
#  Cases

w = np.array([1, 2, 3, 4, 5])
v = np.matmul(A(0.35),w) #<<<<< exact (w, v) couple
v = np.array([1, 2.3, 3.4, 4.5, 5.6])
v = np.array([1, 2.003, 3.004, 4.005, 5.006])


def obj(a):
    Av = np.matmul(A(a),v)
    return np.dot(w,w) - 2*np.dot(w,np.matmul(A(a),v)) + np.dot(Av,Av)

def dobj(a):
    Av = np.matmul(A(a),v)
    dAv = np.matmul(dA(a),v)
    return -2*np.dot(w,np.matmul(dA(a),v)) + 2*np.dot(Av,dAv)    


n = 200
mu = 0.0001
print([obj(a0), a0])
X = [a0]
Y = [obj(a0)]

for j in range(n):
    a1 = a0 - mu*dobj(a0)
    X.append(a1)
    Y.append(obj(a1))
    print([obj(a1), a1])
    a0 = a1

plt.plot(X,Y)
plt.show()