2
$\begingroup$

As written in the title, I have a function $f: \mathbb{R}^n \rightarrow \mathbb{R}$ and I want to find the direction of least change. However, I encountered an inconsistency between two different methods which I will outline below.

I start with the directional derivative as given by this answer: $(\nabla f)\cdot v$. This will give the rate of change along any direction $v = \sum_i\alpha_i\vec{e}_i$ with $\sum_i\alpha_i^2=1$. To obtain the coefficients of the vector corresponding to the least change, I want to minimize the square of that expression: $$ \min_{\{\alpha_i\}} \left[(\nabla f)\cdot v\right]^2 $$

This is similar to the following quadratic problem: $$ \min_{v} \quad v^T \cdot \underbrace{ \begin{pmatrix} \partial_1\partial_1 & \partial_1\partial_2 & \cdots & \partial_1\partial_n \\ \partial_2\partial_1 & \partial_2\partial_2 & \cdots & \partial_2\partial_n \\ \vdots & \vdots & \ddots & \vdots \\ \partial_n\partial_1 & \partial_n\partial_2 & \cdots & \partial_n\partial_n \\ \end{pmatrix} }_{\equiv A} \cdot v \qquad\mathrm{subject}\;\mathrm{to}\quad \lVert v \rVert = 1 $$ where $\partial_i$ is the $i$-th partial derivative evaluated at some reference point.

In order to solve this, I am following the paper William W. Hager, "Minimizing a Quadratic Over a Sphere" (full text available here). Relevant are Lemma 2.1 and 2.2. In summary, the result is that the direction of least change is given by the eigenvector which corresponds to the smallest eigenvalue of $A$. However, it doesn't claim that the smallest eigenvalue must be positive and, in fact, Lemma 2.1 states that only $A - \lambda_1 I$ be positive-semidefinite ($\lambda_1$ being the smallest eigenvalue and $I$ the identity matrix).

The above findings are in agreement with this answer. That answer links to a Wikipedia article about the Structure Tensor, which seems to be similar to the above matrix $A$ when the window function $w(r)$ is the Dirac delta function. This article, however, states that all eigenvalues of $A$ should be positive (see here; this is also mentioned in this paper, paragraph following eq.(16)). Also the 3D interpretation with an ellipsoid seems to require $\lambda_i\geq 0$ as they are associated with the semi-axes of the ellipsoid. This would imply that $A$ is already positive-semidefinite which seems to conflict with the statement in the above paper.

When I apply the above method, i.e. diagonalizing $A$, to a specific problem ($n=36$), then I find indeed that some of the eigenvalues are negative; $A - \lambda_1 I$ on the other hand is positive-semidefinite. So I'm wondering whether this method is correct or if there is something I'm missing? Does it have to do with the dimensionality of the problem (in my specific case $n=36$)?

The second thing I'm wondering about is that the rate of change is zero if $v$ is chosen perpendicular to $\nabla f$. For $n\geq 3$ there exist infinitely many vectors $v$ that are perpendicular to any other given vector, so this seems to contradict the above finding that the direction of least change is precisely given by one of the eigenvectors.


Numerical example

The following are the 36 elements of the gradient, separated by commas, for my specific problem (it's a complex physics simulation, hence I'm just including the gradient here):

-4.629630012686902774e+01,-1.625366508051229175e+01,-3.641531520770513453e+00,3.641531520770513453e+01,1.354472090042690979e+01,2.953193245502916398e+00,2.629008122312370688e+01,1.239008895481674699e+01,5.195843755245732609e+00,-2.098321516541545861e+01,-8.260059303211164661e+00,-1.776356839400250465e+00,1.050270981295398087e+01,1.243449787580175325e+00,-6.439293542825907934e-01,4.269917752708352054e+01,1.769695501252499525e+01,5.773159728050814010e+00,-2.420286193682841258e+00,1.931788062847772380e+00,1.509903313490212895e+00,-1.521005543736464460e+01,-7.527312106958561344e+00,-2.908784324517910136e+00,3.606004383982508443e+01,1.207922650792170316e+01,2.686739719592878828e+00,2.680078381445127889e+01,1.294520046712932526e+01,4.773959005888173124e+00,-2.051692149507289287e+01,-8.038014698286133353e+00,-1.731947918415244203e+00,9.325873406851314940e+00,1.620925615952728549e+00,-5.995204332975845318e-01

I use SymPy, specifically Matrix.diagonalize, in order to obtained the eigenvalues and eigenvectors:

import numpy as np
from sympy import Matrix

gradient = np.array([-4.629630012686902774e+01,-1.625366508051229175e+01,-3.641531520770513453e+00,3.641531520770513453e+01,1.354472090042690979e+01,2.953193245502916398e+00,2.629008122312370688e+01,1.239008895481674699e+01,5.195843755245732609e+00,-2.098321516541545861e+01,-8.260059303211164661e+00,-1.776356839400250465e+00,1.050270981295398087e+01,1.243449787580175325e+00,-6.439293542825907934e-01,4.269917752708352054e+01,1.769695501252499525e+01,5.773159728050814010e+00,-2.420286193682841258e+00,1.931788062847772380e+00,1.509903313490212895e+00,-1.521005543736464460e+01,-7.527312106958561344e+00,-2.908784324517910136e+00,3.606004383982508443e+01,1.207922650792170316e+01,2.686739719592878828e+00,2.680078381445127889e+01,1.294520046712932526e+01,4.773959005888173124e+00,-2.051692149507289287e+01,-8.038014698286133353e+00,-1.731947918415244203e+00,9.325873406851314940e+00,1.620925615952728549e+00,-5.995204332975845318e-01])
A = Matrix(np.multiply.outer(gradient, gradient))
P, D = A.diagonalize(reals_only=True, sort=True, normalize=True)

print('----- Eigenvalues -----')
print(np.diag(D), end='\n\n')

which gives the following output:

----- Eigenvalues -----
[-2.12609948489760e-13 -1.52827730730385e-13 -1.00673222408694e-13
 -7.83376902220698e-14 -5.00864285330806e-14 -4.16036346697355e-14
 -3.07234097627031e-14 -2.42707863964688e-14 -2.00785824570343e-14
 -1.15863543324252e-14 -1.00286213791179e-14 -4.01745072568617e-15
 -1.65034450852433e-15 -8.84471783248045e-16 -7.34601713232325e-16
 -3.49931683094853e-16 -1.87968744502435e-16 -3.39521826408545e-17
 2.39349396454803e-78 9.63838283004546e-17 4.83310171616047e-16
 5.70723316062587e-16 1.21622134384129e-15 1.67127258435284e-15
 2.52716079698100e-15 3.56914743743705e-15 6.33591850166460e-15
 9.44400674233662e-15 2.04608434272058e-14 2.93396570040894e-14
 4.82240672003627e-14 5.21039117277875e-14 7.82644482744230e-14
 1.11237991654078e-13 2.17480897644168e-13 10853.3563960944]

The first 35 eigenvalues have very small magnitude (compared to the 36th one), some of them are negative and others positive.

$\endgroup$
9
  • 1
    $\begingroup$ Rather than looking at the square, you could simply note that $\nabla f \cdot v = 0$ is a linear system of equations on the entries of $v$. $\endgroup$ Jan 20, 2022 at 16:55
  • $\begingroup$ @BenGrossmann Right, that's the other inconsistency which I forgot to include in my original question. $\nabla f\cdot v = 0$ implies that any vector $v$ which is perpendicular to the gradient is a direction of zero change, i.e. the minimum change possible. For $n \geq 3$ there are infinitely many ways to choose such a vector which seems to contradict the other solution where the direction of least change is given by one of the eigenvectors (i.e. one specific direction). $\endgroup$
    – a_guest
    Jan 20, 2022 at 19:11
  • $\begingroup$ Ultimately, the question comes down to what you mean by "the rate of change along a direction". If all you're looking for is a direction $v$ such that we have $$ \frac{d}{dt} f(x + tv) = \lim_{t \to 0} \frac{f(x + tv) - f(x)}{t}= 0, $$ this holds iff $v \cdot \nabla f(x) = 0$. Hager's approach seems to consider a "second order" rate of change. That is, I believe that the goal in that context is to find a $v$ such that $$ \lim_{t \to 0}\frac{f(x + tv) - f(x)}{t^2} $$ is minimized. Note that in that paper, "minimum" does not refer to minimum absolute value. $\endgroup$ Jan 20, 2022 at 22:36
  • $\begingroup$ Normally, this "second order rate" is only considered when $\nabla f = 0$. $\endgroup$ Jan 20, 2022 at 22:39
  • $\begingroup$ Actually, the paper that you linked has no explicit connection to calculus; the goal seems to be merely to minimize $x^TAx - 2b^Tx$ without any notion of looking for a "direction of least change". $\endgroup$ Jan 20, 2022 at 22:44

1 Answer 1

3
+50
$\begingroup$

Your matrix $A$ is an outer product of the gradient vector, and therefore a rank-1 PSD matrix. As given in this answer - https://math.stackexchange.com/a/55166/443030 - $A$ has an eigenvalue $0$ of multiplicity $n-1$, and only one eigenvalue equal to $\lVert \nabla f(x) \rVert^2$. The negative values associated with your matrix are in fact numerical inaccuracies, and can be assumed zero given the magnitude of the maximum eigenvalue.

As stated in the paper, the minimizer $v$ would be an eigenvector corresponding to the smallest eigenvalue of $A$, which in this case is $0$. Note, also, that the eigenvector corresponding to the only nonzero eigenvalue is $\nabla f(x)$, which is the direction of maximum change (in fact, any vector that is scalar multiple of the gradient points in the direction of most change). So, this is consistent to the fact that we would expect the least change to occur in the direction of all vectors in the nullspace of the gradient, which is spanned by all the eigenvectors corresponding to the eigenvalue $0$.

$\endgroup$
1
  • $\begingroup$ Thanks, then it makes all sense, since the paper specifies that the solution is any linear combination with norm 1 of all eigenvectors that correspond to the smallest eigenvalue which also spans a $(n-1)$-dimensional subspace. $\endgroup$
    – a_guest
    Jan 25, 2022 at 21:08

You must log in to answer this question.

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