5
$\begingroup$

I have a symmetric positive definite matrix (hessian) $H$ which is unknown and expensive to compute explicitly (circa 30*30)

Indirectly in my code I have a growing list of pairs of unit vectors $u_i$, and noisy estimates of the matrix product $\hat{v_i}=Hu_i + \epsilon_i$, where $\epsilon_i$ can be assumed to be normal i.i.d variables. Note that due to the noise component, occasionally $u_i^TH\hat{v_i}<0$

I was wondering if there is a method to estimate $H$ from the $(u_i,\hat{v_i})$ pairs, which forces $H$ to be S.P.D, with a starting guess $\hat{H_0} = c*I$.

Obviously due to the small size of the vectors, I can store hundreds, potentially thousands, of these vector pairs as the algorithm runs. Potentially earlier pairs might be less relevant as the algorithm continues, so I may only use the last $n$ pairs for any estimation, or attempt to recency weight them (somehow).

As a point of reference, it's part of an optimization process, and so have tried quasi methods like BFGS, however the noise component of $\hat{v_i}$ simply causes it to blow up, or estimate poorly.

Apologies if this is a well known topic, I have tried searching the site and googling but I'm not even sure the term I should be searching for!

Improve this question
$\endgroup$
1
  • $\begingroup$ See my answer below. As for your larger optimization problem in which BFGS blows up due to errors in gradient differences, if bootstrapping can be applied to the gradient samples, then you can bootstrap the BFGS update to choose a BFGS update bootstrap sample which does not blow up, per my manuscript "Honey, I Shrunk the Stochastic Quasi-Newton Update by Selecting My Favorite Bootstrap Gradient Difference Sample", March 2023, lnkd.in/ed2_YSQF $\endgroup$
    – Mark L. Stone
    Commented Jun 25, 2023 at 15:47

1 Answer 1

Reset to default
1
$\begingroup$

Formulate and solve this as a Linear SDP (a.k.a. Linear Matrix Inequality) constrained linear least squares problem. This is a convex, easy to solve (30 by 30 is not very big) optimization problem, for which there are a number of solvers available, such as Mosek, SeDuMi, SDPT3, and COPT, among others. It can easily be formulated and solved with a convex optimization modeling tool such as CVX, YALMIP, CVXPY, CVXR, or Convex.jl.

min (w.r.t $\hat{H}) $ $\|vec([\hat{v_1} -\hat{H}u_1, ....,\hat{v_n} -\hat{H}u_n])\|_2 $

s.t. $\hat{H} - \text{m}I \succeq 0$

where m is the desired minimum eigenvalue for $H$.

The data points can be weighted if old observations should count less.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you, this is hugely helpful & exactly what I was looking for. I've made some excellent progress with this. So far I'm using CVXPY, and I'm finding convergence to be reliable, but the constraints to be rather 'rubbery'. I.e. H-10I>0 might occasionally throw eigenvalues of 3 or 4 instead of 10+. Have you found one of the packages to be superior to the others in terms of convergence? Additionally if I add, for instance, a ridge regression type penalization to the problem, eg u0,1,2,3,..n=unit vectors in each direction, v0,1,2,...n=lambda*unit direction, it adds instability $\endgroup$
    – SRB121
    Commented Jun 26, 2023 at 11:21
  • $\begingroup$ Did you declare H to be symmetric? Did CVXPY/solver claim to have solved to optimality? If so, you should not have gotten such large infeasibility on eigenvalue unless numerics (problem scaling) are horrible, such as huge input data values, or your syntax was incorrect and you didn't actually specify the problem you intended. $\endgroup$
    – Mark L. Stone
    Commented Jun 26, 2023 at 14:37

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .