Quadproj quickstart

This page gives you a quick overview of the use of quadproj.

The following code snippets should be run sequentially (e.g., to inherit from import), the whole script is available in test_README.

The basics

A simple N dimensional example

In the following snippet, we create an object quadproj.quadrics.Quadric. A quadric is a quadratic hypersurface of the form

\[\mathbf{x}^t A \mathbf{x} + \mathbf{b}^t \mathbf{x} +c = 0.\]

It is created by providing a dict containing the matrix A, the vector b and the scalar c. After definition of the latter parameters, the object is created in a single line (line 17).

We then create some point x0 that we project onto the quadric.

 1from quadproj import quadrics
 2from quadproj.project import project
 3
 4
 5import numpy as np
 6
 7
 8# creating random data
 9dim = 42
10_A = np.random.rand(dim, dim)
11A = _A + _A.T  # make sure that A is positive definite
12b = np.random.rand(dim)
13c = -1.42
14
15
16param = {'A': A, 'b': b, 'c': c}
17Q = quadrics.Quadric(param)
18
19x0 = np.random.rand(dim)
20x_project = project(Q, x0)
21assert Q.is_feasible(x_project), 'The projection is incorrect!'

For data of dimension larger than three, it is unfortunately rather difficult to plot the quadric.

Visualise the solution

Let us try a simpler 2D case.

 1from quadproj.project import plot_x0_x_project
 2import quadproj.utils
 3from os.path import join
 4
 5import pathlib
 6
 7
 8root_folder = pathlib.Path(__file__).resolve().parent.parent
 9output_folder = join(root_folder, 'output')
10
11import matplotlib.pyplot as plt
12
13show=False
14
15A = np.array([[1, 0.1], [0.1, 2]])
16b = np.zeros(2)
17c = -1
18Q = quadrics.Quadric({'A': A, 'b': b, 'c': c})
19
20x0 = np.array([2, 1])
21x_project = project(Q, x0)
22
23fig, ax = Q.plot(show=show)
24plot_x0_x_project(ax, Q, x0, x_project)
25plt.savefig(join(output_folder, 'ellipse_no_circle.png'))

The output is given in output/ellipse_no_circle.png and the projection of x0 onto the quadric (x_project) is depicted as a red point. Mathematically we write the projection \(\mathrm{P}_\mathcal{Q}(x^0)\).

_images/ellipse_no_circle.png

Remark that you can obviously change the output to the current directory (or whatever path).

>>> plt.savefig('ellipse_no_circle.png')

Looking at the previous image, one may believe that the returned projection (red point) is actually not the closest point to x0. This false impression is due to the unequal axes. As a way to verify that the projection is indeed the closest point, we run the following snippet.

1fig, ax = Q.plot()
2plot_x0_x_project(ax, Q, x0, x_project, flag_circle=True)
3fig.savefig(join(output_folder, 'ellipse_circle.png'))
4if show:
5    plt.show()

Switching the parameter flag_circle to True, we verify the optimality of x_project.

_images/ellipse_circle.png

Indeed, we observe that the circle centered in x0 with radius \(\| \mathbf{x} - \mathbf{x}^0 \|\) is:

  • tangent at x_project = \(\| \mathrm{P}_{\mathcal{Q}} (\mathbf{x}^0) \|_2\) ;

  • does not cross the quadric.

Quid of multiple solutions?

The projection is not unique in general! We denote as degenerate case the cases where x0 is located on one of the principal axes of the quadric.

In this degenerate case, it is possible (though not always the case) that multiple points are optima.

For constructing a degenerate case we can:

  • Either construct a quadric in standard form, i.e., with a diagonal matrix A, a nul vector b, and c=-1 and define x0 with a least one entry equal to zero;

  • Or choose any quadric and select x0 to be on any axis.

Let us illustrate the second option. We create x0 by applying the (inverse) standardization from some x0 with at least one entry equal to zero.

Here, we chose to be close to the center and on the longest axis of the ellipse so as to be sure that there are multiple (two) solutions.

Note that the program returns only one solution. Multiple solutions is planned in future releases.

1x0 = Q.to_non_standardized(np.array([0, 0.1]))
2x_project = project(Q, x0)
3fig, ax = Q.plot(show_principal_axes=True)
4plot_x0_x_project(ax, Q, x0, x_project, flag_circle=True)
5fig.savefig(join(output_folder, 'ellipse_degenerated.png'))
6if show:
7    plt.show()

The output figure ellipse_degenerated.png is given below. It can be shown that the reflection of the x_project along the largest ellipse axis (visible because show_principal_axes=True) yields another optimal solution.

_images/ellipse_degenerated.png

Supported quadrics

Ellipse

See previous section for examples of ellipses.

Hyperbola

We illustrate a degenerated projection onto an hyperbola.

In this case, there is no root to the nonlinear function (graphically, the second axis does not intersect the hyperbola). This is not an issue because two solutions are obtained from the other set of KKT points (How does quadproj works?).

1A[0, 0] = -2
2Q = quadrics.Quadric({'A': A, 'b': b, 'c': c})
3x0 = Q.to_non_standardized(np.array([0, 0.1]))
4x_project = project(Q, x0)
5fig, ax = Q.plot(show_principal_axes=True)
6plot_x0_x_project(ax, Q, x0, x_project, flag_circle=True)
7fig.savefig(join(output_folder, 'hyperbola_degenerated.png'))
8if show:
9    plt.show()
_images/hyperbola_degenerated.png

Ellipsoid

Similarly as the 2D case, we can plot ellipsoid.

_images/ellipsoid.png
 1dim = 3
 2A = np.eye(dim)
 3A[0, 0] = 2
 4A[1, 1] = 0.5
 5
 6b = np.zeros(dim)
 7c = -1
 8param = {'A': A, 'b': b, 'c': c}
 9Q = quadrics.Quadric(param)
10
11
12fig, ax = Q.plot()
13
14fig.savefig(join(output_folder, 'ellipsoid.png'))
15
16Q.get_turning_gif(step=4, gif_path=join(output_folder, Q.type+'.gif'))

To ease visualisation, the function get_turning_gif lets you write a rotating gif such as:

_images/ellipsoid.gif

One-sheet hyperboloid

You can easily create a turning gif of a modification of the plots (e.g., by adding the projected point) by having a look at the core of get_turning_gif.

 1
 2
 3A[0, 0] = -4
 4
 5param = {'A': A, 'b': b, 'c': c}
 6Q = quadrics.Quadric(param)
 7
 8x0 = np.array([0.1, 0.42, -1.5])
 9
10x_project = project(Q, x0)
11
12fig, ax = Q.plot()
13ax.grid(False)
14ax.axis('off')
15plot_x0_x_project(ax, Q, x0, x_project)
16ax.get_legend().remove()
17
18save_gif = True
19if save_gif:
20    quadrics.get_gif(fig, ax, elev=15, gif_path=join(output_folder, Q.type+'.gif'))
21if show:
22    plt.show()
_images/one_sheet_hyperboloid.gif

But wait! It seems that the red point is not the closest isn’t it?!

Fortunately, this is not a bug but rather the same scaling issue from Visualise the solution and plotting the iso-curve of the distance function (flag_circle=True), we see that the red ellipsoid (i.e., a distorted ball) is tangent to the quadric at x_project.

1fig, ax = Q.plot()
2ax.grid(False)
3ax.axis('off')
4plot_x0_x_project(ax, Q, x0, x_project, flag_circle=True)
5ax.get_legend().remove()
6
7if save_gif: 
8    quadrics.get_gif(fig, ax, elev=15, gif_path=join(output_folder, Q.type+'_circle.gif'))
_images/one_sheet_hyperboloid_circle.gif

Two-sheet hyperboloid

A Two-sheet hyperboloid is a 3D quadratic surface with two positive eigenvalues and one negative.

 1A[0, 0] = 4
 2A[1, 1] = -2
 3A[2, 2] = -1
 4
 5b = np.array([0.5, 1, -0.25])
 6
 7c = -1
 8
 9param = {'A': A, 'b': b, 'c': c}
10Q = quadrics.Quadric(param)
11
12x0 = np.array([0.1, 0.42, -0.45])
13
14x_project = project(Q, x0)
15
16fig, ax = Q.plot()
17ax.grid(False)
18ax.axis('off')
19plot_x0_x_project(ax, Q, x0, x_project)
20
21if save_gif:
22    quadrics.get_gif(fig, ax, gif_path=join(output_folder, Q.type+'.gif'))
23if show:
24    plt.show()
_images/two_sheet_hyperboloid.gif