Is my Maximize expression too complicated for Mathematica to handle?

I am new to Mathematica and am using it to determine when a function reaches a maximum.

Maximize[
{x^2 j + y^2 k + z^2 l + x y m + x z n + y z o, x^2 + y^2 + z^2 == 1},
{x, y, z}]

That is, given some constants j, k, l, m, n, and o, solve for x, y, andz` in the unit sphere for where the expression reaches a global maximum.

I initially ran this on the online Development Platform, but it fails with the error "This computation has exceeded the time limit for your plan". I then purchased a Developer plan which raised the computation time limit to 20 minutes, and while it doesn't print an error when attempting to solve, it simply doesn't print anything (so, I assume it also reaches the time limit). I then started a trial of Mathematica desktop thinking that if I could just let it run on my local computer, surely it would eventually complete. After five or six days, I had to restart my computer because of an update.

My question is, is this function too complicated for even a program such as Mathematica to solve? Perhaps the function can be rewritten or a command other than Maximize used to achieve the same result? I know enough calculus to do simple integrals and derivatives and such, but can't recognize if something is not practically solvable by a machine.

Not that it matters for the purpose of solving it, but here is a long-winded explanation for how I arrived at the function and why I am interested in solving it.

EDIT:

The linked document explains what the constants are and how they are calculated, but it is a bit detailed so I'll summarize it here.

My program has a list of three dimensional points that are unit vectors - i.e. they all lie on the surface of a unit sphere positioned at (0, 0, 0). To efficiently find the nearest N points among the list closest to another arbitrary point, the program creates a data structure that sorts the list of points into a hierarchy. Each level of the hierarchy defines a plane that halves the unit sphere, the idea being that points in the list below the plane are on one side of the descendant hierarchy, and points above are on the other side. This type of data structure is called a binary space partition, or BSP. As with all BSPs, determining an optimal halving plane is a non-trivial process, hence my attempt to determine if a mathematical function can be used.

The function itself was created by combining two heuristics, one to emphasize that points are evenly spread on either side of the halving plane, and another to emphasize that the spread is large. These are both determined by using the dot product of each point and the unit normal of the plane (given by x, y, and z in the above function).

As such, solving the function for x, y, and z will give the program a good (best?) halving plane at one level of the hierarchy, meaning the plane with the normal [x, y, z] will split the points efficiently. I don't know if it will mean much to know exactly how each constant is calculated, but if P is the list of points and N is the number of points in the list, some preliminary values are calculated as such:

$$a = \displaystyle\sum_{i=1}^N P_{ix}$$

$$b = \displaystyle\sum_{i=1}^N P_{iy}$$

$$c = \displaystyle\sum_{i=1}^N P_{iz}$$

$$d = \displaystyle\sum_{i=1}^N P_{ix}^2$$

$$e = \displaystyle\sum_{i=1}^N P_{iy}^2$$

$$f = \displaystyle\sum_{i=1}^N P_{iz}^2$$

$$g = 2 * \displaystyle\sum_{i=1}^N P_{ix}P_{iy}$$

$$h = 2 * \displaystyle\sum_{i=1}^N P_{ix}P_{iz}$$

$$i = 2 * \displaystyle\sum_{i=1}^N P_{iy}P_{iz}$$

Then the constants as shown in the function are calculated as:

$$j = d - a^2$$

$$k = e - b^2$$

$$l = f - c^2$$

$$m = g - 2ab$$

$$n = h - 2ac$$

$$o = i - 2bc$$

• What do you know about the six parameters? – Ulrich Neumann Oct 9 '18 at 8:35
• @UlrichNeumann there is no succinct way to answer that, but they are summations (sort of) of x, y, and z components of a list of unit vectors, and the x, y, and z in the function comprise a unit vector that is the normal to a plane that efficiently bisects the list of unit vectors. The idea is that the function, once solved, will be used to optimize creation of a data structure used to query geometric datum. – jorgander Oct 9 '18 at 8:47
• The parameters doesn't depend on x,y,z (" they are summations (sort of) of x, y, and z ")? – Ulrich Neumann Oct 9 '18 at 9:27
• @UlrichNeumann No, I should have been more clear with my use of "x, y, and z". The six parameters are calculated using the components of 3D vectors (the x, y, and z coordinates of each), whereas the x, y, and z symbols in the function are the components of a different unit vector that I am interested in finding. So, the x, y, and z as shown in the function depend on the values of the constants j, k, l, ... , not the other way around. They are not mutually dependent. – jorgander Oct 9 '18 at 17:25

Knowing a bit of linear algebra allows you do to an end-run around the Maximize routine. As noted by @UlrichNeumann, you are trying to maximize the quantity $$\vec{v} \cdot (M \vec{v})$$ over all unit vectors $$\vec{v}$$, where $$M$$ is the matrix $$M = \begin{bmatrix} 2j & m & n \\ m & 2k & o \\ n & o & 2l \end{bmatrix}.$$ This turns out to be equivalent to maximizing the Rayleigh quotient over all vectors. Here's how to do this:

$$M$$ is a real symmetric $$3 \times 3$$ matrix, which means that it has three real eigenvalues $$\lambda_1, \lambda_2, \lambda_3$$ and three corresponding orthonormal real eigenvectors $$\hat{e}_1, \hat{e}_2, \hat{e}_3$$. For each such eigenvector, we have $$M \hat{e}_i = \lambda_i \hat{e}_i.$$ Since the vectors $$\{ \hat{e}_1, \hat{e}_2, \hat{e}_3 \}$$ form a basis for the vector space, any vector $$\vec{v}$$ must be expressible as $$\vec{v} = a_1 \hat{e}_1 + a_2 \hat{e}_2 + a_3 \hat{e}_3$$ with $$a_1^2 + a_2^2 + a_3^2 = 1$$ (since $$\vec{v} \cdot \vec{v} = 1$$.) This means that $$\vec{v} \cdot (M \vec{v}) = \lambda_1 a_1^2 + \lambda_2 a_2^2 + \lambda_3 a_3^2.$$

In particular, if $$\lambda_3$$ is the largest of the eigenvalues, then we can rewrite the above as $$\vec{v} \cdot (M \vec{v}) = (\lambda_1 - \lambda_3) a_1^2 + (\lambda_2 - \lambda_3) a_2^2 + \lambda_3,$$ and since both $$(\lambda_1 - \lambda_3)$$ and $$(\lambda_2 - \lambda_3)$$ are negative, it is fairly evident that this is maximized when $$a_1 = a_2 = 0$$. Thus, the quantity $$\vec{v} \cdot (M \vec{v})$$ is maximized when $$\vec{v}$$ is the normalized eigenvector of $$M$$ with the highest eigenvalue.

Mathematica can, of course, find these eigenvectors and eigenvalues:

M={{2 j, m, n}, {m, 2 k, o}, {n, o, 2 l}};
evalues = Eigenvalues[M]
evectors = Eigenvectors[M]

If you run this on symbolic inputs, you'll get a horrendously long output in terms of Root objects. However, if you assign numerical values to all six variables in this expression, it will return a set of eigenvalues and eigenvectors quickly. You'll want to examine the output of Eigenvalues to find the largest eigenvalue; the corresponding entry of Eigenvectors will then contain the corresponding eigenvector, which will be the value of $$\vec{v}$$ that maximizes $$\vec{v} \cdot (M \vec{v})$$. This is easy enough to do "by hand", but I believe the following automated code will also do the trick:

Take[evecs, First[Position[evals, Max[evals]]]]

In conclusion, everyone who wants to do linear programming should know a bit of linear algebra.

• You are right, I don't know enough linear algebra to disseminate all of this, but I don't doubt that your answer is correct, or at least will help me reach a solution. I hope you don't mind that I accept @UlrichNeumann reply since as you say his answer is correct and he provided it first, but many thanks to you for explaining it in more detail. – jorgander Oct 9 '18 at 19:59
• @Michael Seifert: thanks for your clear explanations which reminds me on topics liek Raleigh quotient & Co. – Ulrich Neumann Oct 10 '18 at 7:07
• @UlrichNeumann: I'd never heard this called the "Rayleigh quotient", but it makes sense that the technique would be named. I've inserted a link for further reading. – Michael Seifert Oct 10 '18 at 12:49
• @ MichaelSeifert Thanks for your reply! – Ulrich Neumann Oct 10 '18 at 12:58

The functional

J = ((x^2)*j) + ((y^2)*k) + ((z^2)*l) + (x*y*m) + (x*z*n) + (y*z*o);

you want to maximize is a quadratic form

M={{2 j, m, n}, {m, 2 k, o}, {n, o, 2 l}};
Simplify[1/2 {x, y, z}.M.{x, y, z} == J]
(* True*)

Therefore the maximum depends on the properties of the matrix M. But M depends on the parameters! Without further knowledge of the unknown parameters it seems to be hopeless to find a general solution.

Reduce gives a bunch of conditions:

Reduce[{M.{x, y, z} == 0, x^2 + y^2 + z^2 == 1}, {x, y, z}]
• Thanks for your answer but I'm not sure if my reply to your earlier comment addresses it or not. Specifically, the values of the constants j, k, l, ... do not depend on x, y, and z. – jorgander Oct 9 '18 at 18:37