Fun with shear operations and SVD – II – Shearing of rectangles and cubes with Python and Matplotlib

Posted on 18. August 2023 by eremo

This post series is about shear transformations and related matrices.

In the end we want to better understand the effect of shear operations on “multivariate standard normal distributions” [MSND]. A characteristic property of a MNSD is that its contour hyper-surfaces of constant density are nested spheres. We want to understand the impact of a shear operation on such a random vector distribution.

In the first post of this series
Fun with shear operations and SVD – I – shear matrices and examples created with Blender
we have already had a preliminary look on shear transformations. We have seen that shear transformations are a special class of affine transformations. They are represented by square unipotent, upper (or lower) triangular matrices. The eigenspace is 1-dimensional and the only eigenvalue is 1.

We have also understood that a shear operation can not be reduced to a simple scaling operation with different factors by choosing some special rotated Euclidean coordinate system [ECS], whose origin coincides with volume and symmetry center of the original figure.

Due to its matrix properties which are very different from orthogonal matrices a shear operation is neither a rotation, nor a scaling. Instead we suspect that it might only be decomposed into a combination of a rotation plus a scaling – whatever coordinate system we choose. We still have to prove this in a forthcoming post.

In this post we want to explicitly apply shear matrices to simple 2D- and 3D-figures, namely squares and cubes. We use Python to control the operations and visualize the results with the help of Matplotlib. We will check that the position and orientation of all extremal points in y-direction (2D-case) and z-direction (3D-case) are preserved during shearing. We will depict this by applying the transformations to position-vectors, which point to elements of distinct vertical layers. In 3D I also test the effects of shearing on cubes by just using the 8 corner points of the cube to define 6 faces, which then get transformed.

Afterward we briefly look at the chaining of shear operations and the resulting product of their respective matrices. As some people may think that one could use combinations of shears induced by a mix of upper triagonal and lower triagonal matrices, we investigate the effects of such combinations. We see that such combinations hurt conditions we have already defined for pure shear operations. Even a symmetric matrix, which has the same shear factors in the upper and lower triangular parts of the matrix has a different effect than a real shearing: It distorts a body in the sense of a pure scaling with different factors in potentially all orthogonal directions.

Hint: I omit Python code in this post. If I find the time I will deliver some code snippets at the end of this post series. However, it is really easy to set up your own programs to reproduce the results of my experiments and respective plots with the help of Numpy and Matplotlib.

Before showing results of numerical experiments I first list up some general affine properties and other specific properties of shear transformations.

Affine properties of shear transformations

Important properties of all affine transformations are:

Collinearity: points on a line are transformed into points on a line.
Parallelism: Parallel lines remain parallel lines. Lines remain lines. Planes are transformed into planes.
Preservation of sub-dimensionality: Affine transformations preserve the dimensionality of objects confined to a sub-space of the “affine” space they operate on.
Convexity: Convex figures remain convex – i.e. the sign of a curvature radius does not change.
Proportionality: Distance ratios on parallel line segments are preserved.

Extremal points on convex surfaces may remain extremal with respect to distances to the origin, but not with respect to a specific coordinate axis (rotation!). Note that the third point holds for convex (hyper-) surfaces – as e.g. that of a sphere or ellipsoid.

Properties of shear transformations and their matrices M_sh

I list some properties which I have already discussed in my previous post:

Fix point: When we omit a translation it is natural to choose an Euclidean coordinate system [ECS] such, that its origin coincides with the symmetry center of the body to be sheared. This point then remains a fixed point of the transformation.
M_sh is a unipotent matrix with non-zero-elements only in its upper (or lower) triangular part. The matrix elements on the diagonal are just 1.
det(M_sh) = 1. A shear transformation is invertible. Its rank is rank(M_sh) = n.
The only eigenvalue ε is ε = 1. Its multiplicity in the ℝⁿ is n.
The eigenspace ES_sh has a dimension dim(ES_sh) = 1.
A single pure shear mediated by a unipotent upper (or lower) matrix can not be reduced to a pure scaling operation in some cleverly chosen, rotated coordinate system.

For reasons of convenience we below again chose an Euclidean coordinate system [ECS] whose origin coincides with the center of volume or the symmetry center of a body which we apply a shear transformation to. This not only guarantees a central fixed point. The structure of the diagonal matrix and the 1s on the diagonal

\[ \begin{align} \pmb{\operatorname{M}}_{sh} \, &= \,
\begin{pmatrix}
1 & m_{12} &\cdots & m_{1n} \\
0 & 1 &\cdots & m_{2n} \\
\vdots &\vdots &\ddots &\vdots \\
0 & 0 &\cdots & 1 \end{pmatrix}, \\
\pmb{v}^{sh} \: &= \: \pmb{\operatorname{M}}_{sh} \circ \pmb{v} \phantom{\Huge{(}}
\end{align}
\]

make it clear that the last component of any affected vector remains constant during a shear transformation:

\[
\begin{align}
\pmb{v} \: &= \: (\, v_1, \, v_2, \,\cdots \, v_n)^T, \\
\pmb{v}^{sh} \: &= \: (\, v_1^{sh}, \, v_2^{sh}, \,\cdots \, v_n^{sh})^T \phantom{\huge{(}} \\
v^{sh}_n \: &= \: 1.0 \, *\, v_n \phantom{\huge{(}}
\end{align}
\]

This is independent of the linear coupling between other components. Therefore, an important feature of a shear operation is:

Extremal points on the surface of a sheared body in the n-th coordinate direction remain extremal and are moved within a hyper-plane orthogonal to the sub-space ℝ¹ of the ℝⁿ.

And, more generally:

Points of a sheared body lying in a (n-1)-dimensional sub-space of the ℝⁿ, which are defined by vⁿ = const, are moved within the sub-space, only, during shearing.

In 3D with coordinates (x, y, z) such a sub-space is a hyper-plane defined by z = const.. E.g. for a cube the points of the top face in z-direction are elements who just move within this plane during a shearing-operation. The same holds for any layer of points defined by other values z = const. Actually, the extremal points in z-direction even move along defined straight lines within the plane. In the general case these lines cross the (n-1)-dimensional sub-space.

It is time to take Python to perform some numerical experiments in 2D and 3D – and visualize the above statements. As said, in this post we concentrate on the shearing of squares and cubes. And the behavior of extremal points and other points located on layers with a constant value of the preserved vector component. I admit that squares and cubes are a bit boring, but this post is just a preparational step towards a closer investigation of the the impact of shear operations on statistical vector distributions controlled by Gaussians.

2D: Shearing of a square/rectangle

The figure below shows the end-points of position-vectors, which define data points on a rectangular meshes within a square. All the points reside inside and on the boundary of the square and reflect horizontal layers of the square:

Now we apply a shear matrix with the following elements [[1.0, 1.5], [0.0, 1.0]]. As we only work in two dimensions our matrix just has 4 elements with values ≠ 0. We name the undisturbed vectors of our square v_sq. The resulting vectors after the shear operation will be called v_sh. I.e.:

\[ \pmb{\operatorname{M}}_{sh} \: = \: \begin{pmatrix} 1.0 & 1.5 \\ 0.0 & 1.0\end{pmatrix}
\]

and

\[ \pmb{v}_{sh} \: = \: \pmb{\operatorname{M}}_{sh} \circ \pmb{v}_{sq}
\]

A plot of all the v_sh gives us:

Exactly, what we expected. We directly see the constant shift of layers against each other. And we also see the linearity of the overall effect in the diagonal lines between the new edges. The overall result is a parallelogram.

The points on all the layers were moved along hyper-surfaces with y = const., which in our 2D case are just lines parallel to the x-axis.

Note that the symmetry center at the origin of our figure remained fixed. Point symmetries with respect to the original symmetry center have been preserved. However, original symmetries with respect to the coordinate axes have been broken.

For a rectangle we get analogous, stretched images:

Shear transformation of a cube

The last plots were simple, but visualized the very nature of a shear transformation: Linear differences of (position) vectors between layers. Breaking of plane symmetries for polytopes, while point symmetries are preserved. Let us try to demonstrate something similar in 3D. Once again an illustrative method to do this is to produce multiple planes on different z-levels filled with points on a regular mesh. During plotting we allow for some transparency of these points.

All images show the situation before shearing. You see the flat squares (extending in x- and y-direction) on 4 different z-levels. It can be seen that the dimensions are -15 ≤ v_i ≤ 15, for i ∈ {1, 2, 3}. We now apply a shear matrix M_sh = [[1.0, 0.7, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]

\[ \pmb{\operatorname{M}}_{sh} \: = \: \begin{pmatrix} 1.0 & 0.7 & 0.0
\\ 0.0 & 1.0 & 0.0 \\
0.0 & 0.0 & 1.0
\end{pmatrix}
\]

to the position-vectors of our layer points. I.e. we only induce a linear effect of the y-components onto the x-components of our vectors. The result is:

We get a systematic y-dependent change of the x-coordinate values of our points. The z-level of our layer points, however, remain unchanged. Now, we repeat our experiment, but this time we add an additional shear via a linear dependency of x- and y-component values of our vectors on their z-component by

\[ \pmb{\operatorname{M}}_{sh} \: = \: \begin{pmatrix} 1.0 & 0.7 & 1.2
\\ 0.0 & 1.0 & 1.0 \\
0.0 & 0.0 & 1.0
\end{pmatrix}
\]

The shift of the layers against each other in x- and y-direction gets a clearly visible z-dependence:

But, again and of course, all elements of our 4 layers remain on their respective z-level.

Shearing of a cube: Surface layers

An alternative method to visualize the impact of a shear-transformation on a cube with the help of Matplotlib is based on 8 vectors that span the 6 surface elements of the cube. We also make the cube-faces a bit transparent in some plots. The original cube then looks like:

The result for our last shear matrix

\[ \pmb{\operatorname{M}}_{sh} \: = \: \begin{pmatrix} 1.0 & 0.7 & 1.2
\\ 0.0 & 1.0 & 1.0 \\
0.0 & 0.0 & 1.0
\end{pmatrix}
\]

is:

Here wee see again that the z-level of the upper and lower limiting plates of the cube remained constant. They were just moved in x- and y-direction:

Our shear operation has produced a parallelepiped or a spat. This is consistent to the results of Blender experiments depicted in my previous post.

Remarks on a significant deficit of Matplotlib for 3D

A significant problem that one stumbles into with plots like those above is that Matplotlib does not offer us a real 3D engine. Meaning: It often calculates the so called z-order of objects like hyper-planes (or parts of such), curved manifolds, bodies along our fictitious line of visibility wrongly. So, the shading of e.g. transparent and opaque surfaces according to overlaps in 3D along your line of sight is almost always depicted in a wrong way.

For complex bodies this can drive you crazy and makes you wonder if it would be better to learn controlling Blender with Python. For our relatively simple cases I solved the z-order problem by disabling the computed z-order

ax = fig.add_subplot(111, projection=’3d’, computed_zorder=False)

and by manually assigning a z-order-value to each of the cube’s faces.

Still: As soon as we distort a simple body as a cube strongly (see below) you need a really good overview about the order of your surfaces and their colors after rotations. It is worthwhile to make sketches or to create a similar body with Blender first to clearly see, what happens during a sequence of rotations.

Combinations of shear transformations

Shearing transformations can, of course, be chained in an orderly sequence of operations. Let us take a 2D example of two shear matrices with the same shear-parameter:

\[ \begin{align} \pmb{\operatorname{M}}_{sh} \circ \pmb{\operatorname{M}}_{sh} \: &= \:
\begin{pmatrix} 1 & \lambda \\ 0 & 1 \end{pmatrix} \circ
\begin{pmatrix} 1 & \lambda \\ 0 & 1 \end{pmatrix}
\\
&= \: \begin{pmatrix} 1 & 2\lambda \\ 0 & 1 \end{pmatrix} \phantom{\Huge{(}^T}
\end{align}
\]

For different shear parameters we get

\[ \begin{align} \pmb{\operatorname{M1}}_{sh} \circ \pmb{\operatorname{M2}}_{sh} \: &= \:
\begin{pmatrix} 1 & \lambda \\ 0 & 1 \end{pmatrix} \circ
\begin{pmatrix} 1 & \mu \\ 0 & 1 \end{pmatrix}
\\
&= \: \begin{pmatrix} 1 & \lambda \, +\, \mu \\ 0 & 1 \end{pmatrix} \phantom{\Huge{(}^T}
\end{align}
\]

I.e., we just increase the shearing effect by performing two shear operations after another. If we had taken the negative value (μ = – λ) we would just have inverted the shearing. Note:

A chaining of pure shear matrices, which we had defined as upper triangular matrices in the first post, is a symmetric operation with respect to the order of the matrices.

Combination of upper and lower triangular shear matrices?

However, what happens if we combine an unipotent upper triangular shear matrix with a lower triangular matrix with the same shearing coefficient?

\[ \begin{align} \pmb{\operatorname{M}}_{sh} \circ \pmb{\operatorname{M}}_{sh}^T \: &= \:
\begin{pmatrix} 1 & \lambda \\ 0 & 1 \end{pmatrix} \circ
\begin{pmatrix} 1 & 0 \\ \lambda & 1 \end{pmatrix}
\\
&= \: \begin{pmatrix} 1 \,+\, \lambda^2 & \lambda \\ \lambda & 1 \end{pmatrix} \phantom{\Huge{(}^T}
\end{align}
\]

We get a symmetric matrix. This is not surprising, as this always is the case for product of a matrix with its transposed counterpart. But, in addition, a scaling is superimposed. Note also that the combined operation is not symmetric regarding the order of matrix application:

\[ \begin{align} \pmb{\operatorname{M}}_{sh}^T \circ \pmb{\operatorname{M}}_{sh} \: &= \:
\begin{pmatrix} 1 & 0 \\ \lambda & 1 \end{pmatrix} \circ
\begin{pmatrix} 1 & \lambda \\ 0 & 1 \end{pmatrix}
\\
&= \: \begin{pmatrix} 1 & \lambda \\ \lambda & 1 \,+\, \lambda^2 \end{pmatrix} \phantom{\Huge{(}^T}
\end{align}
\]

Let us apply the second variant to our square (see above) for λ = 1.5:

As expected we get a stronger distortion in y-direction. Note also that in contrast to a normal shear operation z-levels of our 4 level-rows were not preserved. Instead the individual levels appear rotated and stretched. Actually, in one of the next posts we will see that our combined matrix can be factorized into a combination of a rotation, then pure scaling, then back-rotation (by the same angle used in the first rotation). But this corresponds to nothing else than a pure scaling in a rotated ECS.

The funny thing with chaining is that using a simple symmetric matrix from the start for λ = 1.5

\[ \pmb{\operatorname{M}}_{symm} \: = \:
\begin{pmatrix} 1 & \lambda \\ \lambda & 1 \end{pmatrix}
\]

would give us a different result than a sequence of M_sh^T * M_sh:

Note:

Combining a shear matrix with its transposed is something different than just using a matrix with filling in (transposed) elements in the lower part with the same value as their mirrored elements in the upper triangular part.

By combining upper and lower triangular shear matrices we leave the group of pure shear operations

We could call the application of a symmetric matrix a kind of symmetric shearing. But even with a symmetric matrix we loose the central property of a normal shear operation, namely that the z-levels of body-layers parallel to the x/y-plane are preserved. All in all we get into trouble as soon as we combine original shear matrices and their transposed counterparts. Regarding my definition of a shear matrix in my previous post, we obviously leave the group of shear operations by such a combination. Some reader may regard this as an exaggeration. But it makes sense:
Even the symmetric matrix of the combined operation can be reduced to a pure scaling operation in a rotated coordinate system. And this simply is not shearing.

Just for fun: “Symmetric shearing” applied to a cube

Just for the fun of it, let us apply a symmetric version of a 3D-shearing matrix to our cube:
shear = [[1.0, 0.6, 1.0], [0.6, 1.0, 1.2], [1.0, 1.2, 1.0]]

\[ \pmb{\operatorname{M}}_{symm} \: = \: \begin{pmatrix}
1.0 & 0.6 & 1.0 \\
0.6 & 1.0 & 1.2 \\
1.0 & 1.2 & 1.0
\end{pmatrix}
\]

We expect all layers to be stretched into all coordinate directions now; z-levels of the original layers parallel to the x/y-plane will not be preserved. The normal vector to the transformed layers will be inclined to all axes after the operation.

The following plots were done with intransparent, opaque surfaces. Due to the extreme stretching we would otherwise loose the overview – which is nevertheless difficult. I have kept the order of colors for the cubes side-layers the same as for the unstretched cube: sf_color=[‘yellow’, ‘blue’, ‘red’, ‘green’, ‘magenta’, ‘cyan’]. I first give you some selected perspectives:

The distortion in all directions is obvious. Note that the point-symmetries of all corners with respect to the origin are preserved. The following sequence reflects views from a circle around the vertical axis:

What we learn is that even scaling alone can distort a symmetric body quite a lot.

Conclusion

In this post I have visualized the effect of shear operations on squares and cubes. We saw that working with points on layers and with respective vectors is helpful when we want to demonstrate a central property of real shear operations: One vector component stays as it is. Nevertheless, the symmetry breaking of shearing transformations was clearly visible. We also saw that we should represent shear transformations either by upper or lower triangular matrices, but we should not mix these types of matrices.

As a side effect we have seen that a combination of shear matrix with its transposed induces strong deformations in specific coordinate directions by enhanced scaling factors. Symmetric matrices also lead to relativ strong distortions. We have claimed that these distortions can be explained by pure scaling operations in a rotated coordinate system. A proof has to be delivered yet.

In the next post
Fun with shear operations and SVD – III – Shearing of circles
we take a first step in the direction of shearing multivariate normal vector distributions. We keep it simple and just study the impact of shearing transformations on circles.

Fun with shear operations and SVD – I – shear matrices and examples created with Blender

Posted on 14. August 2023 by eremo

In Machine Learning samples of data vectors often undergo linear transformations. Typical examples are the adaption of vector component values due to the choice of a different shifted and/or rotated coordinate system. Another example is the scaling of vectors. As we have seen in an investigation of Autoencoders they produce multivariate normal distributions which actually are the result of linear transformations of elementary Gaussian distributions. And let us not forget that the transport between layers of neural networks includes linear transformations.

In most cases we can interpret a linear transformations in the ℝⁿ in a geometrical way. Linear transformations are members of a sub-class of affine transformations. Affine transformations in turn can be regarded as a combination of

a translation,
a scaling in either coordinate direction,
a reflection = mirroring operation,
a rotation
and/or a shear operation.

Translations are trivial. We also understand the effects of rotation and scaling transformations intuitively. A mirror operation can be understand as a kind of inverted scaling with factors -1.

Scaling operations with different scaling factors along different coordinate directions destroy certain rotational symmetries of a 3-dimensional or n-dimensional body. But scaling transformations do not always eliminate symmetries with respect to mirroring planes: If and when we first choose an Euclidean coordinate system whose axes are aligned with symmetry axes of a body with plane symmetries, we are able to keep up at least some of the planar symmetries of the body during scaling operation.

Something that is a bit more difficult to grasp is a shear operation. The reason is that a shear operation destroys one or multiple plane symmetries of an impacted figure – whatever coordinate system we choose. And a shear operation cannot be reduced to a pure scaling operation in some cleverly chosen coordinate system. Neither can it be reduced to a rotation.

This alone is a good reason to have a closer look at this special kind of affine transformation. Another reason is to study what impact a shear operations has on the construction or transformation of a multivariate normal distribution – a topic which is at the center of another post series in this blog.

In this post I first want to discuss what a matrix that controls a shear operation looks like. We discuss some of its mathematical properties. Afterward I want to demonstrate what effects shearing has on two simple objects – a cube and a sphere. These examples will be done with Blender and the available mouse driven shear operations there.

Blender performs the required math in the background – and the examples below only provide visual impressions.

To get a better understanding of the nature of shear operation we, in the second post of this series, switch to Python and apply shear operations via explicit matrix operations onto vectors. We start with 2D-figures. Afterward we turn to 3D-cuboids and find that we produce spats. Another important example (not only with respect to multivariate normal distributions) is then given by the shearing of 3D-spheres: We will find that such an operation always produce ellipsoids. We will check this mathematically.

In a fourth post I want to discuss an important result from linear algebra, namely the so called SVD-decomposition of a linear transformation and its meaning in our context of shear operations. SVD tells us that we can always replace a shear operation by a sequence of simpler operations. We will test this statement for our 2D- and 3D-objects objects. We then also find that the creation of an ellipsoid at any orientation can be done by applying just one shear operation to a unit sphere. Alternatively, we can perform a scaling of a sphere and rotate it once. This will help us to better understand the creation of multivariate normal distributions in a parallel post series.

Shear operations break symmetries of the affected figures

A shear operation M_sh introduces an asymmetry into originally symmetric figures. It does so by coupling vector components in a specific way. For reasons of simplicity we reduce our argumentation to the ℝⁿ. We describe objects in this space with the help of position vectors to points inside and on the surface of these objects. We specify the components of position vectors with respect to an Euclidean coordinate system [ECS].

To make things easier we assume that the origin of our ECS coincides with the symmetry center of the body we apply the shear operation to. If the body has plane and point symmetries then the origin is placed at the cutting lines or the cutting point of such planes; i.e., the ECS origin coincides with the symmetry-center of the body and coordinate planes will coincide with some of the symmetry planes. For the cases we look at – spheres, cubes, coboids – the ECS origin thus gets identical with the body’s volume center.

A shear transformation applied to a vector couples at least two of the vector’s components such that a growing value for the coordinate value in the direction of a chosen coordinate axis (e.g. v_z) impacts the coordinate value in another coordinate direction (e.g. v_x) in a linear way.

Here is an example for the x- and z-components of a 3-dimensional vector (v_x, v_y, v_z)^T:

\[ \begin{align}
v_x^{sh} \, &= \, v_x \, + \, \alpha_{xz} * v_z \\ v_y^{sh} \,&=\, v_y \\ v_z^{sh}\,&=\, v_z \end{align}
\]

\[ \pmb{\operatorname{M}}_{sh}: \quad \pmb{v} \,=\, \begin{pmatrix} v_x \\ v_y \\ v_z \end{pmatrix} \:\Rightarrow \:
\pmb{v}_{sh} \,=\, \begin{pmatrix} v_x^{sh} \, = \, v_x \, + \, \alpha_{xz} * v_z \\ v_y^{sh} \,=\, v_y \\ v_z^{sh}\,=\, v_z \end{pmatrix},
\]

With α being some constant. This operation breaks an original plane symmetry of a body with respect to the x-z-plane.

A shear operation preserves the orientation of tangent planes at some surface points

But note also, that a shear operation has another unique property that keeps up a geometrical property of the affected body:

A shear operation preserves the orientation of tangent planes to the body’s surface at some particular points. The most important of these points define extrema on the surface with respect to the component whose value is linearly added in the shearing transformation. The tangent planes there do not only keep up their orientation, but also their position.

In the case defined above we would look a extrema with respect to the z-component:

\[ z \,=\, f(x,\,y), \quad { \partial f \over \partial x } \,=\, { \partial f \over \partial y } \, = \, 0
\]

Actually, we would need a bit of multivariate calculus to derive this statement for general bodies. But the example of a cube helps to grasp the central point: The transformation defined above would keep up the orientation and position of the top and the bottom square surfaces of the cube. See below for respective graphics.

In the case of a general n-dimensional body whose surface is closed, continuous and differentiable at all surface points we will find at least two such points whose tangent planes keep their orientation and position during a shearing operation.

Coupling of more than two vector components – and a restriction

Of course we could have introduced more than just one linear coupling of a selected vector component to other components. In three or n dimensions we can couple 2 components to a third one. In n dimensions we can add up shearing contributions of (n-1)-coordinate values to each of the components. But we introduce a restriction:

To avoid confusion and double shearing effects we always set α_zx = 0. And analogously for the other components:

\[ \begin{align} \alpha_{ij} \, &= \, const. \, \neq \, 0, \quad \forall\, i > j \\ \alpha_{ij} \, &= \, 0, \quad \forall\, i < j \end{align} \]

Once again: We avoid a kind of double coupling between the same pair of components. E.g., v_x may become impacted by v_z, but v_z then shall not get impacted by v_x at the same time.

Obviously, this is a an asymmetric policy in the formal handling of vector component coupling. But, it has geometrical consequences (see below). Of course, we follow this asymmetric policy also for other pairs of vector components and respective coupling parameters. This gives us a number of n*(n – 1)/2 potential constant shearing parameters.

A first example produced with the help of Blender

To get a visual impression of what shearing does to a body with some symmetries, let us look at a simple example: a cube.

Regarding the first image we define that the x-coordinate direction goes from the left to the right, the y-direction from the front to the back and the z-direction from the bottom vertically upwards. The last image from the top shows that the x- and y-side-length of our cube are equal. Among other symmetries, a cube and also cuboids obviously have a mirror symmetry with respect to planes parallel to the x/y-plane, the y/z-plane and the x/z-plane through their centers of volume.

Now let us apply two shear operations which couple the x/z-coordinates and the y/z-coordinates. Then we get something like the following:

The first two image show the sheared object from the front, but with the camera moved to different heights. We see that the cube was sheared in two directions – in the x- and in the y-direction. The z-values of the upper and lower surfaces remain constant during the operation. As a result the object gets inclined in two directions: We get a spat (or parallelepiped) as the result of applying a shear operation to a cube or cuboid.

The last image shows the object from exactly the same top position as before the shearing. We see that the original symmetry is destroyed with respect to all 3 orthogonal planes (parallel to coordinate planes) through the center of volume of our object.

Also all of the following objects are the results of simple shearing operations. While this clear for the parallelepiped (derived from a cuboid) and the inclined cylinder (derived from a straight cylinder), it is not directly obvious for the ellipsoid which in the depicted case actually is the result of a shearing applied to a sphere (with x being impacted by y and z). In all cases rotations and translations have been applied in addition to arrange the scene.

We get the impression that some simple geometric shapes like a spat or an ellipsoid could be the direct result of one or more shearing operations applied to originally simple objects with multiple symmetries – as cubes and spheres.

Shear transformations are mathematically done by unipotent upper triangular matrices

Let us return to a more mathematical description of a shear operation. As any other linear operation applied to vectors of the ℝⁿ we can represent a shear operation by a matrix. We just have to express the linear coupling between two vector components by a suitable arrangement of our shear-constants as matrix-elements. By experimenting a bit we find that a square matrix which only has non-zero elements on its diagonal and in its upper right triangular part does the job. We speak of an “upper triangular matrix“:

\[ \pmb{\operatorname{M}} \, = \,
\begin{pmatrix}
m_{11} & m_{12} &\cdots & m_{1n} \\
0 & m_{22} &\cdots & m_{2n} \\
\vdots &\vdots &\ddots &\vdots \\
0 & 0 &\cdots & m_{33} \end{pmatrix}
\]

That such a matrix mixes vector components in the intended way follows directly from the definition of standard matrix multiplications. (Note: In Python/Numpy we talk about “matmul”-operations.)

Note that not all elements in the upper triangular part (aside the diagonal) are required to be different from zero. But at least one should be ≠ 0.

Not a full and neither a symmetric matrix

Now we come back to a point introduced already above: In principle we could allow for a coupling of two selected vector components in both logically possible ways. E.g., when we have an impact of the z- onto the x-component, why not also allow for an analogous impact of x- onto z at the same time? Maybe with a different linear factor?

Well, we then would get a full linear operation with possibly all n**2 elements having different values or we would get a symmetric matrices. A full generalization is not our objective. We want a special sub-case of a linear transformation. But why not allow for a symmetric matrix?

The reason for this will become clear when we discuss the decomposition and factorization of shear matrices in future posts. For the time being the following hint may help: The complexity of a linear operation can be measured by the minimum number of elementary operations required in a well chosen ECS. Can we find a coordinate system where the transformation reduces to just one elementary operation (translation, scaling, reflection, rotation)?

For a shearing transformation this would have to be a scaling operation as this kind of operation can at least potentially break rotational and planar symmetries.

Now we remember a result of linear algebra: Any symmetric matrix corresponds to an operation for wich we always can find a rotated ECS (at the center of an originally highly symmetric figure) in which the transformation of the figure reduces to a pure scaling operation.

For a shear operation we even want to exclude this kind of reduction: We shall not be able to find a ECS in which the operation just becomes a scaling. (A shearing shall in any coordinate system at least require a scaling and a rotation.) How can we achieve this? Well, we must not even allow for a symmetric matrix. For bodies which show original symmetries with respect to planes through its center this means that a shear operation destroys at least some of its planar symmetries whatever coordinate system we choose. Despite the fact that a shear operation can be inverted (see below)!

So, our elimination of some rotational and at least one planar symmetry in an originally highly symmetric body is reflected by the triangular form of the matrix. We avoid a 2-fold coupling between two selected coordinates at the same time. So when we have a matrix element m_i,j, with j > i, we always set the corresponding element to zero: m_j,i = 0.

This restricts the distortions we are going to get. Some surface elements are going to remain congruent and some point symmetries are kept up.

Why did we not use a lower triangular matrix?

Well, this is just a convention. Actually, the transposed matrix of an upper triangular one would also do the job – so we could also take matrices with non-zero elements only in the lower left triangle. But, for convenience, we will stick to matrices with non-zero elements in the upper right triangular part.

What about the diagonal elements of a shear matrix?

First of all note that not all elements on the diagonal of the matrix M_sh should become zero – because then we would eliminate our original figure totally. But for a shear operation we do not want any of the vector components to disappear! In addition: As non-zero elements in the diagonal would reflect a simple scaling we do request that all of the diagonal matrix elements should be equal to 1 for a pure shear operation:

\[ \pmb{\operatorname{M}}_{sh} \, = \,
\begin{pmatrix}
1 & m_{12} &\cdots & m_{1n}(\ne0) \\
0 & 1 &\cdots & m_{2n} \\
\vdots &\vdots &\ddots &\vdots \\
0 & 0 &\cdots & 1 \end{pmatrix}
\]

With at least some of the m_i,j(>i) ≠ 0. We then speak of a unipotent or unitriangular matrix.

An example for 4 dimensions would be:

\[ \pmb{\operatorname{M}}_{sh}^{4D} \: = \: \begin{pmatrix}
1.0 & 0 & 1.5 & 0.8 \\
0 & 1.0 & 0.2 & 0 \\
0 & 0 & 1.0 & 2.0 \\
0 & 0 & 0 & 1.0
\end{pmatrix}
\]

As the diagonal is fixed by elements m_i,i(>i) = 1, we have reduced the number of free and independent (off-diagonal) parameters N^ind_sh(>i) of a pure shear matrix M_sh from n**2 to:

\[ N^{ind}_{sh}(\pmb{\operatorname{M}}_{sh}) :\: = \: n\,*\, (n\,-\,1)/2
\]

So, in 2 dimensions instead of 4 parameters we set the diagonal elements to 1 and choose just 1 free off-diagonal parameter. In 3 dimensions we only have to set 3 instead of 6 off-diagonal parameters.

Basic properties of a shear matrix

A quick calculation shows: The determinant of a shear matrix is directly given by the product of its diagonal elements, and thus we have:

\[ \operatorname{det} \left(\pmb{\operatorname{M}}_{sh}\right) : = \: 1
\]

Therefore, a shear matrix is always invertible. This means its effects can be reverted – also in geometrical terms. The matrix also has full rank r = n.

The eigenspace and the eigenvalue(s) of a shear matrix M_sh are also interesting. We can use a property of triangular matrices to get the eigenvalue(s):

The diagonal elements of a triangular matrix are its eigenvalues.

This in turn can be understood from analyzing the characteristic polynomial det(M – λI) = 0.

Thus: The one and only eigenvalue ev of M_sh is 1.

\[ \pmb{ev}_{\operatorname{M}_{sh}} \, = \, 1
\]

The eigenvalues algebraic multiplicity in the characteristic polynom is n This does not yet tell us much about the geometrical multiplicity, i.e. the dimension of the eigenspace.

But all the n eigenvalues share the same eigenspace ES. ES is the kernel of the matrix (A – λI). Therefore, we can derive the eigenspace’s dimension by the dimension formula for this particular matri:

\[ \operatorname{dim}_{ES} \, = \, n \,- \, \operatorname{rank}\left(\pmb{\operatorname{M}}_{sh} \, – \, \lambda*\pmb{\operatorname{I}} \right) \,=\, n \, – \, (n\,-\,1) \,=\, 1
\]

as the rank of the combined matrix on the right side is (n – 1).

You also see this understand this by looking at the 4-dim example matrix given above:
When you write down the equations you will find that only 1 variable can be chosen freely. So the kernel of the matrix (A – λI) is only 1-dimensional. Logically and geometrically, this is clear as any other coordinate value (≠ 0) beyond the first one mixes in contributions.

It is easy to prove that the eigenspace ES can be based upon the eigenvector e₁.

\[ \pmb{e}_1 \, = \, \left( \, 1,\, 0,\, \cdots,\,0\right)^T, \: \\ \mbox{with}\, (n\,-\, 1) \, \mbox{ zero components}
\]

A shear matrix cannot be diagonalized

Can M_sh be diagonalized? Answer: No! From linear algebra we know that this would require n linearly independent eigenvectors. But we have just one! So:

\[ \pmb{\operatorname{M}}_{sh} : \neq \: \pmb{\operatorname{Q}} \circ \pmb{\operatorname{D}}_{sh} \circ \pmb{\operatorname{Q}^T}
\]

Think a bit about the consequences:
An operation like that on the right side of the un-equation above would correspond to the representation of our shear matrix in a different (rotated) coordinate system. A diagonalization corresponds to an eigendecomposition; in geometrical terms a diagonal matrix represents a pure scaling operation. But, obviously, such a description of our shear operation in some (cleverly) selected ECS is not possible:

We cannot find a (translated and rotated) ECS in which a shear operation reduces to a scaling transformation. In any chosen (rotated) ECS there is always an additional second rotation involved which we cannot simply invert or get rid off!

We shall see this in more detail when we discuss decompositions. We simply cannot get rid of the fundamental asymmetry we have build into M_sh by not allowing non-zero elements in the lower triangular part of the matrix! No change to a whatever selected coordinate system will help.

So, a shear operation really is a special symmetry breaking operation which cannot be reduced to just one elementary operation. In any ECS it requires at least a rotation and a scaling!

Summary:

A pure shear operation is introduced by some off-diagonal and non-zero element in the upper right triangular part of a unipotent matrix. All off-diagonal elements in the lower left triangular part are set to zero.

More illustrations from Blender

First some images for a cuboid that was more extremely deformed by shearing in x- and y-direction than the first one:

Now we come to an interesting example: We apply shearing in x- and y-direction to a sphere. We get something whose shape looks much more regular than that of a cuboid. We pick a unit sphere and apply to shearing operations in x- and y-direction:

The first image shows the original sphere from some point of view. The second image was taken from the same perspective, but now after the shearing. The resulting figure looks very regular – and it seems to have plane symmetries.

Another example with a bit of different lighting gives us:

Actually, we get an ellipsoid. We clearly can identify the two points which have a maximum absolute z-values whose tangent planes did not change tier orientation.

Ellipsoids (in 3D) are figures with 3 main axes and planar symmetries. The original fully rotational symmetry (first image of the series above) is obviously broken. However, we seem to have saved at least some planar symmetries. Hmmm …. Just when we thought we had understood that a shearing operation breaks original planar symmetries we have in our new example case kept some.

Why is this? What is so special about the shearing of spheres? Keep this question in mind.

Conclusion

Shearing transformations are a special class of linear and thus transformations. Shearing obviously can break some plane symmetries of originally highly symmetric bodies – as cubes. In contrast to other elementary affine operations no coordinate system can be found where the number of operations to reproduce the effects of a shearing can be reduced to one. In particular it is not possible to find an Euclidean coordinate system where it is just reduced to a (symmetry breaking) scaling with different factors in different coordinate directions.

Blender offers a simple graphical method to apply shearing to simple geometrical bodies. This allowed us to get an impression of the effects on cubes and spheres. Something that was a bit astonishing was the shearing effect on a unit spheres. It just created ellipsoids – again highly symmetric bodies WITH planar symmetries.

In the next post of this series
Fun with shear operations and SVD – II – Shearing of rectangles and cubes with Python and Matplotlib
we will look closer at shearing operation with the help of Python and matplotlib.

Multivariate Normal Distributions – II – random vectors and their covariance matrix

Posted on 6. August 2023 by eremo

This post series is based on notes I took when I recently studied the mathematical properties of Multivariate Normal Distributions [MNDs]. A motivation to publish the notes came from the analysis of a Machine Learning experiment: A Convolutional Autoencoder [CAE] trained on human face images had produced a MND in its latent space. But MNDs are, of course, relevant in other contexts, too.

In this post we take some preparatory steps: I want to introduce the term “random vector” to describe statistical multivariate distributions in vector form. We then have a look at the definition of a probability density function [pdf] in multiple dimensions. Furthermore we discuss linear transformations of a random vector, its expectation values and its covariance matrix. The latter both marks and quantifies (linear) relations between its components. I will also point out a simple but useful equation for the behavior of the covariance matrix when a matrix operates on its random vector. Eventually we also have a look at marginal distributions and their probability density function.

I will proceed not too formally and my discussion of the topics will not be really complete in a pure mathematical sense. On our way I will also remind you of some basic ideas and equations from statistics. People who are familiar with all the terms named above can jump over this post.

Random Vectors and multivariate distributions

Let us assume that we can describe an object of interest, e.g. in the context of an ML problem, by a set of n indexed, real-value variables s_j. Instead of a manifest object we could also use n variables to characterize the result of a well defined experiment. When I speak of “objects below I include this more abstract option.

Each object can mathematically be represented by a data point in the ℝⁿ. We simply relate each of the indexed coordinate axes of a corresponding Euclidean coordinate system with a certain object property. The coordinates of the data points are then given by the variable values.

Let us now order a complete set of parallel measured property values s_j^c for a concrete object as indexed components of a vector s^c. The vector then represents our object:

\[ \pmb{s}^c \: = \: \left( \, s_1^c, \, s_2^c, \, ….\, s_n^c\, \right)^T \,
\]

The T indicates the transposition operation. Note: Throughout this post series we assume that untransposed vectors are written in vertical form. I.e. the components are arranged as rows.

Samples

Let us assume that we have m concrete observations of individual objects (or experimental results). The values of the n real variables for each object are always measured or qualified in parallel:

We either randomly pick a bunch of m concrete objects out of a wider population of available objects and determine the values of their n characteristic s_j-variables.
Or we repeat a concrete experiment, whose outcome is described by values for a set of n measurable s_j-variables, many times. The setup for our experiment defines precisely what the variables are and how we measure their values. The experimental result can be interpreted as an abstract “object”.

Thus we get samples of object data. Our experimental setup and the properties of our objects define a sample space Ω we have a clear description of how we produce one or multiple vectors s. A particular observation maps Ω to an individual vector.

Let us assume that at least in principle our vector components can take any real value. Producing a vector s can then formally be regarded as a function to the ℝⁿ:

F_S: Ω => ℝⁿ, s = F_S(Ω, experiment number, … ).

A sequence of many observations, i.e. creating a sample, formally corresponds to a sequence of many independent applications of this function. An implicit (and sometime problematic) assumption is that picking multiple objects or our experiments do not influence Ω and the statistics of the basic population.

The outcome, namely a set of resulting s-vectors, depends on the statistical properties of the underlying object population. We index the k-th vector by s^k.

Each sample can be visualized by a bunch of data points in an n-dimensional space. A concrete sample in 3 dimensions may e.g. look like this from two perspectives:

But our sample can also be represented by a data matrix. We can e.g. arrange the columns to be given by the s^k vectors. If we have a sample of m objects then we get a (n x m)- data matrix M_S ∈ ℝ^nxm. A row with index j then provides us with a sample statistics for the j-th object property.

\[
\pmb{\operatorname{D}}_S \: = \: \begin{pmatrix}
s^1_1 & s^2_1 &\cdots & s^m_1 \\
s^1_2 & s^2_2 &\cdots & s^m_2 \\
\vdots &\vdots &\ddots &\vdots \\
s^1_n & s^2_n &\cdots & s^m_n \end{pmatrix}
\]

\[ \pmb{\operatorname{D}}_S \: = \: \left(\, \pmb{s}^1, \pmb{s}^2, \cdots …, \pmb{s}^m\right).
\]

Often you will find the transposed matrix for reasons of an efficient notation. I prefer this version as it later allows for matrix operations on the vector components from the left side, e.g. in a Python program.

Random vectors

We indicate the statistical element in creating a concrete vector from an underpinning population by a vector symbol S. By using this symbol we understand implicitly that any concrete vector S = s^c is created by applying our well defined function F_S.

It is natural to split S up into n individual univariate random variables S_j: Ω => ℝ, which represent the statistics per vector component. We write S in vector form:

\[ \pmb{S} \: = \: \left( \, S_1, \, S_2, \, ….\, S_n\, \right)^T
\]

So:

\[ \pmb{S} \:=\: \pmb{s}^c\: =: \left(\, S_1=s_1, S_2=s_2, \cdots …, S_n=s_n \, \right)^T.
\]

The vector variable S gets concrete component values by a well defined statistical process. But, if you like to see it more globally, S represents any object from the population and thus in a certain meaning the whole population – and its statistics.

We call S a multivariate random variable or a random vector.

Thus, our statistical multivariate distribution can now be represented by a random vector – and we know how to create samples by turning S repeatedly into concrete vectors. Why all the fuss?

Linear operations on random vectors

The vector notation allows us to combine a random vector with mathematical vector operations, in particular linear transformations. We implicitly understand that such an operation is applied to any vector of the population (and in certain contexts also to elements of our samples). Such a transformation consists of a mapping of vectors of the original population or a sample to new vectors in the same or a different vector space by one and the same linear transformation. The mapping gives us new component values.

In linear algebra we learn that a linear vector transformation can be represented by a matrix plus a translation vector. If the transformation maps a random vector S of the ℝ^m to a random vector X into the ℝⁿ then the linear (affine) transformation is defined by a (n x m)-matrix A ∈ ℝ^nxm and a vector b ∈ ℝⁿ:

\[ \begin{align}
\pmb{x} \,&=\, (x_1, x_2, …, x_n)^T \,=\, \pmb{\operatorname{A}} \pmb{s} \, + \, \pmb{b} \\
&\, \\
\pmb{X} \, &= \, \pmb{\operatorname{A}S} + \pmb{b} \end{align}
\]

In close similarity to standard operations in vector spaces a linear transformation of a random vector can have two meanings:

The operation either actively transforms our random vector (and the underlying vector distribution) – e.g. by a rotation in the ℝⁿ.
Or we perform a passive switch to another Euclidean coordinate system (e.g. by a rotation of the coordinate axes) and the transformation adapts the component values accordingly.

Invertible transformations?

If and when we have a transformation within the ℝⁿ the matrix A would typically be a real valued quadratic (n x n)-matrix. In many cases we will assume that the transformation and its matrix A are invertible. This imposes certain conditions on the matrix: A must then be positive definite and have a determinant det(A) ≠ 0.

\[ \operatorname{det}\left(\pmb{\operatorname{A}}\right) \, =\, \left|\pmb{\operatorname{A}}\right|
\, \ne \, 0
\]

But note that in the more general case also rectangular matrices can have a “left”- or a “right”-inverse matrices; see https://en.wikipedia.org/ wiki/ Generalized_ inverse.

Rectangular matrices, e.g. with m (< n) rows and n columns have a rank r ≤ m < n, and not all of the column vectors of the matrix are linearly independent.

Probability density functions in the ℝⁿ

Let us assume that the population behind a random vector covers the ℝⁿ completely and in a continuous way. By “continuous” we mean that a transition to infinitesimal small volumes at any position in the ℝⁿ is possible: We still will find vectors of the population there and component values change smoothly between adjacent elements.

We can apply a frequency analysis to the vector component values of the objects of a huge sample or of a number of samples. E.g. we can define intervals into which the components may fall or not. The analysis gives us a probability P(S=s_ΔV) for the occurrence of position vectors s_ΔV which point into a certain finite region ΔV of the ℝⁿ.

If you better like to think in terms of data points:

We get a probability that a data point representing a statistically picked object lies within a certain defined and finite region of the ℝⁿ.

How we split up the ℝⁿ into separate volume regions depends on available information about the object population or on an analysis of sufficiently large samples.

If each component of the vectors s can take any real value we may under certain circumstances make a transition to probability densities. Probability densities can be distilled from huge samples or alternatively very many samples picked from the underlying population. To get an idea of the density function we, in practice, raise the number of observations and reduce the volume regions of the ℝⁿ to a proper averaging size that avoids wild fluctuations. The latter is called “sampling”.

If we by many experiments can assume a converging limiting process to infinitesimal real value intervals ds_j for each of the vector components then we end up with a bunch of n continuous 1-dim probability density functions p_j(s_j): ℝ => ℝ. A concrete probability value for a component must be derived via integration of p_j() over a finite interval a ≤ s_j ≤ b:

\[ \operatorname{P}(a \le s_j \le b) \,=\, \int_a^b p(s_j)ds_j
\]

An interpretation as a probability requires, of course, a normalization of such an integral over the whole ℝ.

It is clear that the full probability density function p_S(s) of the vector distribution varies continuously with the 1-dim probability density functions of its components.

In this post series we are interested in examples of vector distributions which can be described by a continuous multidimensional probability density function [p.d.f] for infinitesimal volume elements dV_S = ds₁, ds₂…ds_n:

\[ p_{\pmb{S}}(\pmb{s}) \: = \: p_{\pmb{S}} \left( \, s_1, \, s_2, \, ….\, s_n\, \right)
\]

\[ \begin{align} P(\pmb{s}, V) \, &= \, \int_V p_{\pmb{S}}(\pmb{s}) \, ds_1ds_2…ds_n \\
&= \, \int\int…\int_V p_{\pmb{S}} \left( \, s_1, \, s_2, \, ….\, s_n\, \right) \, ds_1ds_2 \cdots ds_n
\end{align}
\]

where V marks some finite volume element in ℝⁿ. Note again: A normalization of such an integral over the whole ℝⁿ is required to get well defined probability values. The normalization factor above is assumed to be integrated in the mathematical expression for p_S(s).

For the purposes of this post series we assume that both p_S and the p_j() are differentiable.

Correlations of random vector components

There may be dependencies which control how the components of a random vector vary together. So the variation of p_S(s) may depend in a rather complicated way on the density functions of their components p_j(s_j)

\[ p_{\pmb{S}}(\pmb{s}) \: = \: p_{\pmb{S}} \left( \, p_1(s_1), \, p_2(s_2), \, ….\, p_n(s_n) \, \right)
\]

The variation of components may e.g. happen in pairwise correlated or uncorrelated ways. I.e. we can in general not assume that the total probability density function p_S(s) for a random vector is e.g. a product of the 1-dimensional density probability functions p_j(s_j for its components.

Such dependencies make the statistics of multivariate vector distributions more complex than the statistics of their constituting (marginal) univariate components. The relation between p_S and the individual p_j can be difficult to derive and it may even reflect non-linear component relations. Luckily we will see that MNDs contain only linear couplings of the components and thus simple correlations described by constant coefficients.

In ML we have only samples with a limited number of individual vectors available. The investigated samples have to be big enough to reveal the properties of the total probability density function p_S(s), which assumedly describes the variation of the underlying continuous vector population in the ℝⁿ.

Probability density after linear transformations and contour hyper-surfaces

A transformation of the components of a random vector has, of course, also an impact on the resulting probability density:

\[ p_{\pmb{X}}(\pmb{x}) \,=\, p_{\pmb{X}} \left( \pmb{\operatorname{A}s} \,+\, \pmb{b} \right) \, = \, p_{\pmb{X}} \left(\, x_1, \, x_2, \, ….\, x_n\, \right)
\]

and for a finite volume V_X in the target space

\[ P\left(\pmb{x}, V_X\right) \,=\, \int\int…\int_{V_X} \, = \, p_{\pmb{X}} \left( x_1, \, x_2, \, ….\, x_n \right)\, dx_1dx_2 \cdots dx_n
\]

In case of a well defined distribution S we can hope to find an analytical expression for the multidimensional probability p_X(x) of the transformed random vector X. But as we will see we have to invest some work into the calculation. And we will answer the question by what parameters the probability density of a MND is defined in a unique way.

We expect that the concrete form of p_X is influenced in a characteristic way by the properties of the original p_j(s_j)-pdfs and, of course, by our transformation matrix A. This is the case for MNDs and we will study the effects of linear transformations on MNDs and their projections to sub-spaces in more detail in forthcoming posts.

As the matrix A mixes components we may also expect the creation of (linear) relations between the X-components. Therefore, the shape of the density distribution may change due to scaling, mirroring, shearing and rotation effects. This results from the fact that linear transformations are affine transformations in geometrical terms.

This leads to the question by what method we can describe the “shape” of a multivariate distribution. A simple method is to determine hyper-surfaces of constant probability density values. I.e. contour-surfaces:

\[ p_{\pmb{X}}(\pmb{x}) \, = \, p_{\pmb{X}} \left(\, x_1, \, x_2, \, ….\, x_n\, \right) \:=\: const.
\]

This actually is the definition of a (n-1)-dimensional hypersurface. The dependency on the component values can be complicated. But for the special case of a MND we will find a rather comprehensible form which we can interpret. So we will get a chance to understand and qualify the effects of a linear transformation by its impact on the contour hyper-surfaces of the probability density of a MND.

Expectation value and covariance of a statistical vector distribution

The covariance of two univariate distributions X and Y is defined as

\[
\begin{align}
\operatorname{cov}(X, Y)
&= \operatorname{\mathbb{E}}\left[\,\left(\,X \,-\, \operatorname{\mathbb{E}}\left(X\right)\,\right) \, \left(\,Y \,-\, \operatorname{\mathbb{E}}\left(Y\right)\, \right)\,\right] \\
&= \operatorname{\mathbb{E}}\left(X Y\right) \,-\, \operatorname{\mathbb{E}}\left(X\right) \operatorname{\mathbb{E}}\left(Y\right)
\end{align}
\]

The expectation value and the covariance of a general vector distribution S are defined by:

\[ \begin{align} \pmb{\mu}\left(\pmb{S}\right) \,=\, \pmb{\mu}_S \: &= \: \operatorname{\mathbb{E}}\left(\pmb{S} \right) \, := \, \left( \operatorname{\mathbb{E}}(S_1), \, \operatorname{\mathbb{E}}(S_2), …, \, \operatorname{\mathbb{E}}(S_n) \right)^T \\
\operatorname{Cov}\left(\pmb{S}\right) \: &:= \: \operatorname{\mathbb{E}}\left[\left(\pmb{S} – \operatorname{\mathbb{E}}(\pmb{S}) \right)\, \left(\pmb{S} – \operatorname{\mathbb{E}}(\pmb{S}) \right)^T \right] \end{align}
\]

Note the order of transposition in the definition of Cov! A (vertical) vector is combined with a transposed (horizontal) vector. The rules of a matrix-multiplication then give you a matrix as the result! And the expectation value has to be taken for every element of the matrix.

Thus, the interpretation of the notation given above is:

Pick all pairwise combinations (S_j, S_k) of the component distributions. Calculate the covariance of the pair cov(S_j, S_k) and put it at the (j,k)-place of the matrix.

Meaning:

\[ \begin{align}
\operatorname{Cov}\left(\pmb{S}\right)\:&=\:\operatorname{\mathbb{E}} {\begin{pmatrix} (S_{1}-\mu_{1})^{2}&(S_{1}-\mu_{1})(S_{2}-\mu_{2})&\cdots &(S_{1}-\mu_{1})(S_{n}-\mu_{n})\\(S_{2}-\mu_{2})(S_{1}-\mu_{1})&(S_{2}-\mu_{2})^{2}&\cdots &(S_{2}-\mu _{2})(S_{n}-\mu_{n})\\\vdots &\vdots &\ddots &\vdots \\(S_{n}-\mu_{n})(S_{1}-\mu_{1})&(S_{n}-\mu_{n})(S_{2}-\mu _{2})&\cdots &(X_{n}-\mu_{n})^{2} \end{pmatrix}}
\\\\
&=\quad {\begin{pmatrix} \operatorname{Var} (S_{1})&\operatorname{Cov} (S_{1},S_{2})&\cdots &\operatorname{Cov} (S_{1},S_{n})\\ \operatorname{Cov} (S_{2},S_{1})&\operatorname{Var} (S_{2})&\cdots &\operatorname{Cov} (S_{2},S_{n})\\\vdots &\vdots &\ddots &\vdots \\ \operatorname{Cov} (S_{n},S_{1})&\operatorname{Cov} (S_{n},S_{2})&\cdots &\operatorname {Var} (S_{n}) \end{pmatrix}}
\end{align}
\]

Keep the details of this definition in mind!

The above matrix is the (variance-) covariance matrix. As Wikipedia tells you, some people call it a bit different. I refer to it later on just as the covariance matrix (of a random vector).

Useful properties of the covariance matrix of a random vector

Note that the diagonal elements of the covariance matrix are just the variances of the individual component distributions, whilst the off diagonal elements contain the pairwise covariances of the vector components. In case of a vector distribution with uncorrelated (or even independent) S_j the matrix becomes diagonal:

\[
\begin{pmatrix}
(\sigma_1^s)^2 & 0 &\cdots & 0 \\
0 & (\sigma_2^s)^2 &\cdots & 0 \\
\vdots &\vdots &\ddots &\vdots \\
0 & 0 &\cdots & (\sigma_n^s)^2 \end{pmatrix}
\]

Other properties can be derived from the linearity of expectation values and the fact that

\[ \begin{align} \operatorname{Cov}\left(aX+bY,\,cW+dV\right) \,=\, ac*\operatorname{Cov}\left(X,\, W\right)\,&+\,ad*\operatorname{Cov}\left(X,\,V\right)\,\\ +\,bc*\operatorname{Cov}\left(Y,\,W\right)\,&+\,bd*\operatorname{Cov}\left(Y,\,V\right) \end{align}
\]

It is relatively easy to show that the defined expectation value of a vector distribution behaves linearly with respect to linear transformations of S with a matrix A:

\[
\operatorname{\mathbb{E}}\left( \pmb{\operatorname{A}S} + \pmb{b} \right) \, = \pmb{\operatorname{A}} \operatorname{\mathbb{E}}\left(\pmb{S} \right) \,+\, \pmb{b} \vphantom{(}
\]

And the covariance matrix Cov transforms as follows:

\[ \begin{align}
\operatorname{Cov}\left(\pmb{\operatorname{A}S}\right) \: &= \: \operatorname{\mathbb{E}}\left[\, \left( \pmb{\operatorname{A}} (\pmb{S} \,-\, \mu_S ) \right) \, \left( \pmb{\operatorname{A}} (\pmb{S} \,-\, \mu_S ) \right)^T \, \right] \\ &= \: \operatorname{\mathbb{E}}\left[\, \pmb{\operatorname{A}} ( \pmb{S} \,-\, \mu_S ) \, (\pmb{S} \,-\, \mu_S)^T \pmb{\operatorname{A}}^T \,\right] \\
&= \: \pmb{\operatorname{A}} \operatorname{\mathbb{E}}\left[\, ( \pmb{S} \,-\, \mu_S ) \, (\pmb{S} \,-\, \mu_S)^T \,\right] \pmb{\operatorname{A}}^T \\
\operatorname{Cov}\left(\pmb{\operatorname{A}S}\right) \: &= \: \pmb{\operatorname{A}} \operatorname{Cov}\left(\pmb{S}\right) \pmb{\operatorname{A}}^T \end{align}
\]

The intermediate steps for the Cov – and in particular the 3rd one – are in my opinion best understood when you just write the stuff down for e.g. (3×3)-matrices and carefully regard the Cov-definition for random vectors. You will see that the definition of Cov just allows for the right index juggling! The extension to more dimensions is straightforward and could formally be done by induction.

See e.g. also https://cds.nyu.edu/ wp-content/ uploads/ 2021/05/ covariance_matrix.pdf and https://stats.stackexchange.com/ questions/ 113700/ covariance-of-a-random-vector-after-a-linear-transformation.

Other useful properties are

\[
\begin{align}
\operatorname{Cov}\left(V_j,\, V_k\right) \,&=\, \operatorname{\mathbb{E}}\left(V_j, \,V_k\right) \,-\, \operatorname{\mathbb{E}}\left(V_j\right)\operatorname{\mathbb{E}}\left(V_k\right) \\
\operatorname{Cov}\left(\pmb{S}\right) \,=\,
\operatorname{Cov}\left(\pmb{S},\, \pmb{S}\right) \,&=\, \operatorname{\mathbb{E}}\left(\pmb{S}, \,\pmb{S}^T\right) \,-\, \operatorname{\mathbb{E}}\left(\pmb{S}\right)\operatorname{\mathbb{E}}\left(\pmb{S}\right)^T
\end{align}
\]

By the way, for two different random vectors X and Y the matrix

\[
\operatorname{Cov}\left(\pmb{X},\, \pmb{Y}\right) \,=\, \operatorname{\mathbb{E}}\left(\pmb{X}, \,\pmb{Y}^T\right) \,-\, \operatorname{\mathbb{E}}\left(\pmb{X}\right)\operatorname{\mathbb{E}}\left(\pmb{Y}\right)^T
\]

is called the cross-covariance matrix.

Remark on independence and correlation

Two univariate random variables X and Y are statistically (!) independent when their joint probability, i.e. their probability for the occurrence of a value pair (x,y) is the product of the individual probabilities for the occurrence of x and y:

\[
p_{X,Y}(x,y) \,=\, p_X(x)*p_Y(y)
\]

However, uncorrelatedness does mean something else, namely that the covariance of the distributions is 0:

\[ \begin{align}
\operatorname{Cov}\left(X,\, Y\right) \,&=\, 0 \\
\Rightarrow \: \operatorname{\mathbb{E}}\left(X Y\right) \,&=\, \operatorname{\mathbb{E}}\left(X\right) \operatorname{\mathbb{E}}\left(Y\right)
\end{align}
\]

Independence and uncorrelatedness are two different (!) things and uncorrelatedness does in general not imply statistical independence!

Thus the uncorrelatedness of a pair of components (S_j, S_k) of a random vector S does unfortunately not automatically imply their independence!

We shall later in this series that this is actually a bit different for a MND.

Conclusion

In this post we have collected and clarified some basic terms which we will use in forthcoming texts about multivariate normal distributions. An important step was to describe a statistical vector distribution in form of a random vector. We also have seen that a random vector can be subject to linear transformations.

We have discussed a limiting transition to probability density functions for statistical vector distributions. We understood that dependencies and correlations may make it difficult to derive the form of the probability density for the random vector from the probability density functions for the vector components.

In addition we have defined the mean value and the covariance matrix for a random vector. We saw how a linear transformation affects the covariance matrix.

In the next post we will apply our knowledge to a simple random vector based on Gaussian probability density functions without correlations. We will study what effects orthogonal and shearing transformations have on such a distribution.

How to plot a tri-variate data distribution in 3D and add surfaces of confidence ellipsoids from a trivariate normal distribution with the same covariance matrix

Posted on 28. July 2023 by eremo

In Machine Learning [ML] or statistics it is interesting to visualize properties of multivariate distributions by projecting them into 2- or 3-dimensional sub-spaces of the datas’ original n-dimensional variable space. The 3-dimensional aspects are not so often used because plotting is more complex and you have to fight with transparency aspects. Nevertheless a 3-dim view on data may sometimes be more instructive than the analysis of 2-dim projections. In this post we care about 3-dim data representations of tri-variate distributions X with matplotlib. And we add ellipsoids from a corresponding tri-variate normal distribution with the same covariance matrix as X.

Multivariate normal distribution and their projection into a 3-dimensional sub-space

A statistical multivariate distribution of data points is described by a so called random vector X in an Euclidean space for the relevant variables which characterize each object of interest. Many data samples in statistics, big data or ML are (in parts) close to a so called multivariate normal distribution [MND]. One reason for this is, by the way, the “central limit theorem”. A multivariate data distribution in the ℝⁿ can be projected orthogonally onto a 3-dimensional sub-space. Depending on the selected axes that span the sub-space you get a tri-variate distribution of data points.

Whilst analyzing a multivariate distribution you may want to visualize for which regions of variable values your projected tri-variate distributions X deviate from adapted and related theoretical tri-variate normal distributions. The relation will be given by relevant elements of the covariance matrix. Such a deviation investigation defines an application “case 1”.

Another application case, “case 2”, is the following: We may want to study a 3-variate MND, a TND, to get a better idea about the behavior of MNDs in general. In particular you may want to learn details about the relation of the TND with its orthogonal projections onto coordinate planes. Such projections give you marginal distributions in sub-spaces of 2 dimensions. The step from analyzing bi-variate to analyzing tri-variate normal distributions quite often helps to get a deeper understanding of MNDs in spaces of higher dimension and their generalized properties.

When we have a given n-dimensional multivariate random vector X (with n > 3) we get 3-dimensional data by applying an orthogonal projection operator P on the vector data. The relatively trivial operator projects the data orthogonally into a sub-volume spanned by three selected axes of the full variable space. For a given random vector their will, of course, exist multiple such projections as there is a whole bunch of 3-dim sub-spaces for a big n. In “case 2”, however, we just create basic vectors of a 3-dim MND via a proper random generation function.

Regarding matplotlib for Python we can use a scatter-plot function to visualize the resulting data points in 3D. Typically, plots of an ideal or approximate tri-variate normal distribution [TND] will show a dense ellipsoidal core, but also a diffuse and only thinly populated outer region. To get a better impression of the spatial distribution of X relative to a TND and the orientation of the latter’s main axes it might be helpful to include ideal contour surfaces of the TND into the plots.

It is well known that the contour surfaces of multivariate normal distributions are surfaces of nested ellipsoids. On first sight it may, hower, appear difficult to combine impressions of such 2-dim hyper-surfaces with a 3-dim scatter plot. In particular: Where from do we get the main axes of the ellipsoids? And how to plot their (hyper-) surfaces?

Objective of this post

The objective of this post is to show that we can derive everything that is required

to plot general tri-variate distributions
plus ellipsoids from corresponding tri-variate normal distributions

from the covariance matrix Σ of our random vector X.

In case 2 we will just have to define such a matrix – and everything else will follow from it. In case 1 you have to first determine the (n x n)-covariance matrix of your random vector and then extract the relevant elements for the (3 x 3)-covariance matrix of the projected distribution out of it.

The result will be plots like the following:

The first plot combines a scatter plot of a TND with ellipsoidal contours. The 2nd plot only shows contours for confidence levels 1 ≤ σ ≤ 5 of our ideal TND. And I did a bit of shading.

How did I get there?

The covariance matrix determines everything

The mathematical object which characterizes the properties of a MND is its covariance matrix (Σ). Note that we can determine a (n x n)-covariance matrix for any X in n dimensions. Numpy provides a function cov() that helps you with this task. The relevant elements of the full covariance matrix for orthogonal projections into a 3-dim sub-space can be extracted (or better: cut out) by applying a suitable projection operator. This is trivial: Just select the elements with the (i, i)- and (i, j)-indices corresponding to the selected axes of the sub-space. The extracted 9 elements will then form the covariance matrix of the projected tri-variate distribution.

In case 2 we just define a (3 x 3)-Σ as the starting point of our work.

Let us assume that we got the essential Σ-matrix from an analysis of our distribution data or that we, in case 2, have created it. How does a (3 x 3)-Σ relate to ellipsoidal surfaces that show the same deformation and relations of the axes’ lengths as a corresponding tri-variate normal distribution?

Creation of a MND from a standardized normal distribution

In general any n-dim MND can be constructed from a standardized multivariate distribution of independent (and consequently uncorrelated) normal distributions along each axis. I.e. from n univariate marginal distributions. Let us call the standardized multivariate distribution SMND and its random vector Z. We use the coordinate system [CS] where the coordinate axes are aligned with the main axes of the SMND as the CS in which we later also will describe our given distribution X. Furthermore the origin of the CS shall be located such that the SMND is centered. I.e. the mean vector μ of the distribution shall coincide with the CS’s origin:

\[ \pmb{\mu} \: = \: \pmb{0}
\]

In this particular CS the (probability) density function f of Z is just a product of Gaussians g_j(z_j) in all dimensions with a mean at the origin and standard deviations σ_j = 1, for all j.

With z = (z₁, z₂, …, z_n) being a position vector of a data point in the distribution, we have:

\[ \begin{align} \pmb{\mu} \: &= \: \pmb{0} \\
f_{\pmb{Z}} \:&=\: \pmb{\mathcal{N}}_n\left(\pmb{0}, \, \pmb{\operatorname{I}}_n \right) \:=\: \pmb{\mathcal{N}}_n\left(\pmb{\mu}=0 ,\, \pmb{\operatorname{\Sigma}}=\pmb{\operatorname{I}}_n \right)
\end{align}
\]

and in more detail

\[ \begin{align} f_{\pmb{Z}}(\pmb{z}) \, =\, \pmb{\mathcal{N}}_n\left(\pmb{0}, \, \pmb{\operatorname{I}}_n \right) \: &= \: \prod\limits_{j=1}^n g(z_j, \, \mu_j=0, \, \sigma_j=1) \\
&= \, {1 \over \sqrt{(2\pi)^n} } \,
{\large e}^{ – \, {\Large 1 \over \Large 2} {\Large \sum\limits_{j=1}^n} \, {\Large z_j}^2 } \\
&=\, {1 \over \sqrt{(2\pi)^n} } \,
{\large e}^{ – \, {\Large 1 \over \Large 2} {\Large \pmb{z}^T\pmb{z}} } \end{align}
\]

The construction recipe for the creation of a general MND X_N from Z is just the application of a (non-singular) linear transformation. I.e. we apply a (n x n)-matrix onto the position vectors of the data points in the ℝⁿ. Let us call this matrix A. I.e. we transform the random vector Z to a new random vector X_N by

\[ \pmb{X}_N \: = \: \pmb{\operatorname{A}} \circ \pmb{Z}
\]

The resulting symmetric (!) covariance matrix Σ_X of X_N is given by

\[ \pmb{\operatorname{\Sigma}}_X \: = \: \pmb{\operatorname{A}} \pmb{\operatorname{A}}^T
\]

Σ determines the shape of the resulting probability distribution completely. We can reconstruct an A’ which produces the same distribution by an eigendecomposition of the matrix Σ_X. A’ afterward appears as a combination of a rotation and a scaling. An eigendecomposition leads in general to

\[ \pmb{\operatorname{\Sigma}}_{X_N} \: = \: \pmb{\operatorname{V}} \pmb{\operatorname{D}} \pmb{\operatorname{V}}^{-1}
\]

with V being an orthogonal matrix and D being a diagonal matrix. In case of a symmetric, positive-definite matrix like our Σ we even can get

\[ \pmb{\operatorname{\Sigma}}_{X_N} \: = \: \pmb{\operatorname{V}} \pmb{\operatorname{D}} \pmb{\operatorname{V}}^T
\]

D contains the eigenvalues of Σ, whereas the columns of V are the components of the eigenvectors of Σ (in the present coordinate system). V represents a rotation and D a scaling.

The required transformation matrix T, which leads from the unrotated and unscaled SMND Z to the MND X_N, can be rewritten as

\[ \pmb{\operatorname{A}}’ \: = \: \pmb{\operatorname{T}} \: = \: \pmb{\operatorname{V}} \pmb{\operatorname{D}}^{1/2}
\]

with

\[ \pmb{\operatorname{\Sigma}}_{X_N} \: = \: \pmb{\operatorname{A}}’ \left(\pmb{\operatorname{A}}’\right)^T
\]

The 1/2 abbreviates the square root of the matrix values (i.e. of the eigenvalues). A relevant condition is that Σ_X is a symmetric and positive-definite matrix. Meaning: The original A itself must not be singular!

This works in n dimensions as well as in only 3.

Creation of a trivariate normal distribution

The creation of a centered tri-variate normal distribution is easy with Python and Numpy: We just can use

np.random.multivariate_normal( mean, Σ, m )

to create m statistical data points of the distribution. Σ must of course be delivered as a (3 x 3)-matrix – and it has to be positive definite. In the following example I have used

\[ \pmb{\operatorname{\Sigma}}_X \:=\: {\begin{pmatrix}
31 & -4 & 5 \\
-4 & 26 & 44 \\
5 & 44 & 85
\end{pmatrix}}
\]

We have to feed this Σ-matrix directly into np.random.multivariate_normal(). The result with 20,000 data points looks like:

We see a dense inner core, clearly having an ellipsoidal shape. While the first image seems to indicate an overall slim ellipsoid, the second image reveals that our TND actually looks more like an extended lens with a diagonal orientation in the CS. =>
Note: When making such §D-plots we have to be careful not to make premature conclusions about the overall shape of the distribution! Projection effects may give us a wrong impression when looking from just one particular position and viewing angle onto the 3-dimensional distribution.

Create some nested transparent ellipsoidal surfaces

Especially in the outer regions the distribution looks a bit fluffy and not well defined. The reason is the density of points drops rapidly beyond a 3-σ-level. One gets the idea that a sequence of nested contour surfaces may be helpful to get a clearer image of the spatial distribution. We know already that contour surfaces of a MND are the surfaces of ellipsoids. How to get these surfaces?

The answer is rather simple: We just have to apply the mathematical recipe from above to specific data points. These data points must be located on the surface of a unit sphere (of the Z-distribution). Let us call an array with three rows for x-, y-, z-values and m columns for the amount of data points on the 3-dimensional unit sphere US. Then we need to perform the following transformation operation:

\[ \pmb{S}_X \: = \: \pmb{\operatorname{T}} \circ \pmb{US}
\]

to get a valid distribution S_X of points on the surface of an ellipsoid with the right orientation and relation of the main axes lengths. The surface of such an ellipsoid defines a contour surface of a TND with the covariance matrix Σ.

After the application of T to these special points we can use matplotlib’s

ax.plot_surface(x, y, z, …)

to get a continuous surface image.

To get multiple nested surfaces with growing diameters of the related ellipsoids we just have to scale (with growing and common integer factors applied to the σ_j-eigen-values).

Note that we must create all the surfaces transparent – otherwise we would not get a view to inner regions and other nested surfaces. See the code snippets below for more details.

SVD decomposition

We also have to get serious about numerically calculating T. I.e., we need a way to perform the “eigendecomposition”. Technically, we can use the so called “Singular Value Decomposition” [SVD] from Numpy’s linalg-module, which for symmetric matrices just becomes the eigendecomposition.

Code snippets

We can now write down some Jupyter cells with Python code to realize our SVD decomposition. I omit any functionality to project your real data into a 3-dim sub-space and to derive the covariance matrix. These are standard procedures in ML or statistics. But the following code snippets will show which libraries you need and the basic steps to create your plots. Instead of a general MND, I will create a TND for scatter plot points:

Code cell 1 – Imports

import numpy as np
import matplotlib as mpl
from matplotlib import pyplot as plt
from matplotlib.colors import ListedColormap
import matplotlib.patches as mpat
from matplotlib.patches import Ellipse
from matplotlib.colors import LightSource

Code cell 2 – Function to create points on a unit sphere

# Function to create points on unit sphere  
def pts_on_unit_sphere(num=200, b_print=True):

    # Create pts on unit sphere     
    u = np.linspace(0, 2 * np.pi, num)
    v = np.linspace(0, np.pi, num)
    
    x = np.outer(np.cos(u), np.sin(v))
    y = np.outer(np.sin(u), np.sin(v))
    z = np.outer(np.ones_like(u), np.cos(v))
    
    # Make array  
    unit_sphere = np.stack((x, y, z), 0).reshape(3, -1)
    
    if b_print:
        print()
        print("Shapes of coordinate arrays : ", x.shape, y.shape, z.shape)
        print("Shape of unit sphere array : ",  unit_sphere.shape)
        print()

    return unit_sphere, x

Code cell 3 – Function to plot TND with ellipsoidal surfaces

# Function to plot a TND with supplied allipsoids based on the sam Sigma-matrix
def plot_TND_with_ellipsoids(elevation=5, azimuth=5,
                             size=12, dpi=96, 
                             lim=30, dist=11, 
                             X_pts_data=[],
                             li_ellipsoids=[], 
                             li_alpha_ell=[], 
                             b_scatter=False, b_surface=True, 
                             b_shaded=False, b_common_rgb=True,
                             b_antialias=True, 
                             pts_size=0.02, pts_alpha=0.8, 
                             uniform_color='b', strid=1,
                             light_azim=65, light_alt=45
                            ):

    # Prepare figure
    plt.rcParams['figure.dpi'] = dpi
    fig = plt.figure(figsize=(size,size))  # Square figure
    ax = fig.add_subplot(111, projection='3d')

    # Check data
    if b_scatter:
        if len(X_pts_data) == 0:
            print("Error: No scatter points available")
            return

    if b_surface:
        if len(li_ellipsoids) == 0:
             print("Error: No ellipsoid points available")
             return
        if len(li_alpha_ell) < len(li_ellipsoids): 
             print("Error: Not enough alpha values for ellipsoida")
             return
        num_ell = len(li_ellipsoids)

    # set limits on axes 
    if b_surface: 
        lim = lim * 1.7 
    else: 
        lim = lim 

    # axes and labels
    ax.set_xlim(-lim, lim)
    ax.set_ylim(-lim, lim)
    ax.set_zlim(-lim, lim)
    ax.set_xlabel('V_1', labelpad=12.0)
    ax.set_ylabel('V_2', labelpad=12.0)
    ax.set_zlabel('V_3', labelpad=6.0)

    # scatter points of the MND
    if b_scatter: 
        x = X_pts_data[:, 0]
        y = X_pts_data[:, 1]
        z = X_pts_data[:, 2]
        ax.scatter(x, y, z, c='r', marker='o', s=pts_size, alpha=pts_alpha) 

    # ellipsoidal surfaces
    if b_surface: 
        # Shaded
        if b_shaded: 
            # Light source 
            ls = LightSource(azdeg=light_azim, altdeg=light_alt)
            cm = plt.cm.PuBu
            li_rgb = []
            for i in range(0, num_ell): 
                zz = li_ellipsoids[i][2,:]
                li_rgb.append(ls.shade(zz, cm))
            
            if b_common_rgb:
                for i in range(0, num_ell-1):
                    li_rgb[i] = li_rgb[num_ell-1]
            
            # plot the ellipsoidal surfaces
            for i in range(0, num_ell): 
                ax.plot_surface(*li_ellipsoids[i], 
                                rstride=strid, cstride=strid, 
                                linewidth=0, antialiased=b_antialias, 
                                facecolors=li_rgb[i], 
                                alpha=li_alpha_ell[i]
                               )
        # just uniform color    
        else:
            print("here")
            #ax.plot_surface(*ellipsoid, rstride=2, cstride=2, color='b', alpha=0.2)
            for i in range(0, num_ell):
                ax.plot_surface(*li_ellipsoids[i], 
                                rstride=strid, cstride=strid, 
                                color=uniform_color, alpha=li_alpha_ell[i])

    ax.view_init(elev=elevation, azim=azimuth)
    ax.dist = dist
    return

Some hints: The scatter data from the distribution X (or X_N) must be provided as a Numpy array. The TND-ellipsoids and the respective alpha values must be provided as Python lists. Also the size and the alpha-values for the scatter points can be controlled. Reducing both can be helpful to get a glimpse also on ellipsoids within the relative dense core of the TND.

Some parameters as “elevation”, “azimuth”, “dist” help to control the viewing perspective. You can switch showing of the scatter data as well as ellipsoidal surfaces and their shading on and off by the Boolean parameters.

There are a lot of parameters to control a primitive kind of shading. You find the relevant information in the online documentation of matplotlib. The “strid” (stride) parameter should be set to 1 or 2. For slower CPUs one can also take higher values. One has to define a light-source and a rgb-value range for the z-values. You are free to use a different color-map instead of PuBu and make that a parameter, too.

Code cell 4 – Covariance matrix, SVD decomposition and derivation of the T-matrix

b_print = True

# Covariance matrix 
cov1 = [[31.0, -4, 5], [-4, 26, 44], [5, 44, 85]]
cov = np.array(cov1)

# SVD Decomposition
U, S, Vt = np.linalg.svd(cov, full_matrices=True)
S_sqrt = np.sqrt(S)

# get the trafo matrix T = U*SQRT(S)
T = U * S_sqrt

# Print some info
if b_print: 
    print("Shape U: ", U.shape, " :: Shape S: ", S.shape)
    print(" S      : ", S)
    print(" S_sqrt : ", S_sqrt)
    print()
    print(" U :\n", U)
    print(" T :\n", T)

Code cell 5 – Transformation of points on a unit sphere

# Transform pts from unit sphere onto surface of ellipsoids 
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
b_print = True
num = 200
unit_sphere, xs = pts_on_unit_sphere(num=num)

# Apply transformation on data of unit sphere
ell_transf = T @ unit_sphere 
# li_fact = [2.5, 3.5, 4.5, 5.5, 6.5]
li_fact = [1.0, 2.0, 3.0, 4.0, 5.0]
num_ell = len(li_fact)
li_ellipsoids = [] 
for i in range(0, num_ell):
    li_ellipsoids.append(li_fact[i] * ell_transf.reshape(3, *xs.shape))

if b_print:
    print("Shape ell_transf  : ", ell_transf.shape)
    print("Shape ell_transf0 : ", li_ellipsoids[0].shape)

Code cell 6 – TND data point creation and plotting

# Create TND-points 
n = 10000
mean=[0,0,0]
pts_data = np.random.multivariate_normal(mean, cov1, n)

# Alpha values for the ellipsoids
common_alpha = 0.05
li_alpha_ell = []
for i in range(0, num_ell):
    li_alpha_ell.append(common_alpha)
#li_alpha_ell = [0.08, 0.05, 0.05, 0.05, 0.05]
#li_alpha_ell = [0.12, 0.1, 0.08, 0.08, 0.08]
li_alpha_ell = [0.3, 0.1, 0.08, 0.08, 0.08]

# Plotting
elev = 5
azim = 15
plot_TND_with_ellipsoids(elevation=elev, azimuth=azim, 
                         lim= 25, dist=11,
                         X_pts_data = pts_data,
                         li_ellipsoids=li_ellipsoids, 
                         li_alpha_ell=li_alpha_ell, 
                         b_scatter=True, b_surface=True, 
                         b_shaded=False, b_common_rgb=False,
                         b_antialias=True, 
                         pts_size=0.018, pts_alpha=0.6,
                         uniform_color='b', strid=2,
                         light_azim=65, light_alt=75
                       )
#

Example views

Let us play a bit around. The ellipsoidal contours are for σ = 1, 2, 3 ,4, 5. First we look at the data distribution from above. The ellipsoidal contours are for σ = 1, 2, 3 ,4, 5. We must reduce the size and alpha of the scatter points to 0.01 and 0.45, respectively, to still get an indication of the innermost ellipsoid of σ = 1. I used 8000 scatter points. The following plots show the same from different side perspectives – with and without shading and adapted alpha-values and number values for the scatter points

Yoo see: How a TND appears depends strongly on the viewing position. From some positions the distribution’s projection of the viewer’s background plane may even appear spherical. But on the other side it is remarkable how well the projection onto a background plane keeps up the ellipses of the projected outermost borders of the contour surfaces (of the ellipsoids for confidence levels). This is a dominant feature of MNDs in general.

A last hint: To get a more volumetric impression it is required to both work with the position of the lightsource, the alpha-values and the colormap. In addition turning off ant-aliasing and setting the stride to something like 5 may be very helpful, too. The next image was done for an ellipsoid with slightly different extensions, only 3 ellipsoids and stride=5:

Conclusion

In this blog I have shown that we can display a tri-variate (normal) distribution, stemming from a construction of a standardized normal distribution or from a projection of some real data distribution, in 3D with the help of Numpy and matplotlib. As soon as we know the covariance matrix of our distribution we can add transparent surfaces of ellipsoids to our plots. These are constructed via a linear transformation of points on a unit sphere. The transformation matrix can be derived from an eigendecomposition of the covariance matrix.

The added ellipsoids help to better understand the shape and orientation of the tri-variate distribution. But plots from different viewing angles are required.

Multivariate Normal Distributions – I – objectives

Posted on 12. July 2023 by eremo

Machine Learning [ML] algorithms are applied to multivariate data: Each individual object of interest (e.g. an image) is characterized by a set of n distinct and quantifiable variables. The variable values may e.g. come from measurements.

A sample of such objects corresponds to a data distribution in a multidimensional space, most often the ℝⁿ. We can visualize our objects as data points in an Euclidean coordinate system of the ℝⁿ: Each axis represents the values a specific variable can take; the position of a data point is given by the variable values.

Equivalently, we can use (position-) vectors to these data points. Thus, when training ML algorithms we typically deal with vector distributions, which by their very nature are multivariate. But also the outputs of some types of neural networks like Autoencoders [AE] form multivariate distributions in the networks’ latent spaces. For today’s ML-scenarios the number of dimensions n can become very big – even if we compress information in latent spaces. For a variety of tasks in generative ML we may need to understand the nature and shape of such distributions.

An elementary kind of a continuous multivariate vector distribution, for which major properties can be derived analytically, is the so called Multivariate Normal Distribution [MND]. MNDs, their marginal and their conditional distributions are of major importance both in the fields of statistics, Big Data and Machine Learning. One reason for this is the “central limit theorem” of statistics (in its vector form).

Some conventional ML-algorithms are even based on the assumption that the population behind the concrete data samples can be approximated by a MND. Due to the central limit theorem we find that averages of big samples of multivariate training data for a population of specific types of observed objects tend to form a MND. But also data samples in latent spaces of neural networks may show a multivariate normal distribution – at least in parts.

For the concrete problem of human face generation via a trained convolutional Autoencoder [CAE] I have actually found that the data produced in the CAE’s latent space can very well be described by a MND. See the posts on Autoencoders in this blog. This alone is motivation enough to dive a bit deeper into the (beautiful) mathematical properties of MNDs.

Just to illustrate it: The following plots show projections of the approximate MND onto coordinate planes.

We find the typical elliptic contour lines which are to be expected for a MND. And here are some generated face images from statistical vectors which I derived from an analysis of the characteristic features of the 2-dim projections of the latent MND which my CAE had produced:

ML is math in the end – and MNDs are no exception

Some of my readers may have noticed that I wanted to start a series on the topic of creating random vectors for a given MND-like vector distribution. The characteristic parameters for the n-dimensional MND can either stem from an analysis of experimental ML data or come from theoretical sources. This was in April. But, I have been silent on this topic for a while.

The reason was that I got caught up in the study of the math of MNDs, of their properties, their marginal distributions and of quadratic forms in multiple dimensions (ellipsoids and ellipses). I had to re-collect a lot of mathematical information which I once (45 years ago) had learned at university. Unfortunately, multivariate analysis (i.e. data analysis in multidimensional spaces) requires some (undergraduate) university math. Regarding MNDs, knowledge both in linear algebra, statistics and vector analysis is required. In particular matrices, their decomposition and their geometrical interpretation play a major role. And when you try to understand a particular problem which obviously is characterized by an overlap of multiple mathematical disciplines the amount of information can quickly grow – without the connections and consistency becoming clear at first sight.

This is in part due to the different fields the authors of papers on MNDs work in and the different focuses they have on properties of MNDs. Although many introductory information about MNDs is available on the Internet, I have so far missed a coherent and comprehensive presentation which illustrates the theoretical insights by both ideal and real world examples. Too often the texts are restricted to pure formal derivations. And none of the texts discussed the problem of vector generation within the limits of MND confidence levels. But this task can become important in generative ML: At high confidence levels outliers become a strong weight – and deviations from an ideal MND may cause disturbances.

One problem with appropriate vector generation for creative ML purposes is that ML experiments deliver (latent) data which are difficult to analyze as they reside in high-dimensional spaces. Even if we already knew that they form a MND in some parts of a latent space we would have to perform a drill down to analytic formulas which describe limiting conditions for the components of the statistical vectors we want to create.

The other problem is that we need a solid understanding of confidence levels for a multidimensional distribution of data points, which we approximate by a MND. And on one’s way to understanding related properties of MNDs you pass a lot of interesting side aspects – e.g. degenerate distributions, matrix decompositions, affine transformations and projections of multidimensional hypersurfaces onto coordinate planes. Far too interesting to refrain from not writing something about it …

After having read many publicly available articles on MNDs and related math I had collected a bunch of notes, formulas and numerical experiments. The idea of a general post series on MNDs grew in parallel. From my own experiences I thought that ML people who are confronted with latent representations of data and find indications of a MND would like to have an introduction which covers the most relevant aspects of MNDs. On a certain mathematical level, and supported by illustrations from a concrete example.

But I will not forget about my original objective, namely the generation of random vectors within confidence levels. In the end we will find two possible approaches: One is based on a particular linear transformation, whose mathematical form is determined by a covariance analysis of our data distribution, and random number generators for multiple Gaussian distributions. The other solution is based on a derivation of precise conditions on random vector components from ellipses which are produced by projections of our real experimental data distribution onto coordinate planes. Such limiting conditions can be given in form of analytic expressions.

The second approach can also be understood as a reconstruction of a multivariate distribution from low-dimensional projection data:

We create vectors of a concrete MND-like vector distribution in n dimensions by only referring to characteristic data of its two-dimensional projections onto coordinate planes.

This is an interesting objective in itself as the access to and the analysis of 2-dimensional (correlated) data may be a much easier endeavour than analyzing the full distribution. But such an approach has to be supported by mathematical arguments.

Objectives of this post series

Objectives of this post series are:

We want to find out what a MND is in mathematical and statistical terms and how it can be based on a simpler vector distributions within the ℝⁿ.
We want to study the basic role of a standardized multivariate normal distribution in the game and the impact of linear affine transformations on such a distribution – in terms of linear algebra and from a geometrical point of view.
We also want to describe and interpret the difference between normal MNDs and so called degenerate MNDs.
We want to understand the most important mathematical properties of MNDs. In particular we want to better grasp the mathematical meaning of correlations between the vector components and their impact on the probability density function. Furthermore the relation of a MND to its marginal distributions in sub-spaces of lower dimensions is of major interest.
We want to formally create a MND-approximation to a real multivariate data distribution by an analysis of real distribution’s properties and in particular from parameters describing the correlations between the vector components. Of particular interest are the covariance matrix and the precision or correlation matrix.
We want to study the role of projections when turning from a MND to its marginal distributions and the impact of such projections on the matrices qualifying the original and its marginal distributions.
We want to understand the form of contour hyper-surfaces for constant probability density values of a MND. We also want to derive what the projections of these hyper-surfaces onto coordinate planes look like.
We want to show that both contour hyper-surfaces of the MND and of its projections in marginal distributions contain the same proportions of integrated data points and, equivalently, the same probability proportions resulting from an integration of the probability density from the distribution’s center up to the hyper-surfaces.
We want to illustrate basic MND-creation principles and the effects of linear affine transformations during the construction process by an ideal 3-dimensional MND example and by projections of a real vector distribution from an ML-experiment onto 2-dimensional and 3-dimensional sub-spaces. We also want to illustrate the relation between the MND and its marginal distributions by plotting concrete 3-dimensional examples and their projections onto coordinate planes.
We want to use the derived MND properties for the creation of statistical vectors v which fulfill the following conditions:
- Each of the generated v is a member of a vector population, which has been derived from a ML experiment and which to a good approximation can be described by a MND (and its extracted basic parameters).
- Each v has an endpoint within the multidimensional volume enclosed by a contour-hypersurface of the MND’s probability density function [p.d.f.],
- The limiting hypersurface is defined by a chosen confidence level.
We want to create statistical vectors within the limit of contour hyper-surfaces by using elementary construction principles of a MND.
In a second approach we want to reduce vector creation to solving a sequence of 2-dimensional problems. I.e. we want to work with 2-dim marginal distributions in 2-dim sub-spaces of the ℝⁿ. We hope that the probability density functions of the relevant distributions can be described analytically and provide computable limiting conditions on vector components.
Note: The production of statistical vectors from data of projected low-dimensional marginal distributions corresponds to a reconstruction of the full MND from its projections.
During random vector creation we want to avoid PCA-transformations of the whole real data distribution or of projections of it.
Based on MND-parameters we want to find analytic expressions for the vector component limits whenever possible.

The attentive reader has noticed that the list above includes an assumption – namely that a multidimensional contour hypersurface of a MND can be associated with something like a confidence level. In addition we have to justify mathematically that the reduction to data of 2-dimensional projections of the full vector distribution is a real option for statistical vector creation.

The last three points are a bit tough: Even if we believe in math textbooks and get limiting hyper-curves of a quadratic form in our coordinate planes the main axes of the respective ellipses may show angles versus the coordinate axes (see the example images above). All this would have to be taken care of in a precise analytic form of the limits which we impose on the components of our aspired statistical vectors.

So, this series is, at least in parts, going to be a tough, but also very satisfactory journey. Eventually, after having clarified diverse properties of MNDs and their marginal distributions in lower dimensional spaces, we will end up with quadratic equations and some simple matrix operations.

Objectives of the next post

We must not forget that statistics plays a major role in our business. In ML we deal with finite collections (samples) of individual object data which are statistically picked from a greater population (with assumed statistical properties). An example is a concrete collection of images of human faces and/or their latent vectors. The data can be organized in form of a two-dimensional data matrix: Its rows may indicate individual objects and its columns properties of these objects (or vice versa). In either direction we have vectors which focus on a particular aspect of the data: Individual objects or the statistics of a specific object property.

While we are used to univariate “random variables” we have to turn to so called “random vectors” to describe multidimensional statistical distributions and respective samples picked from an underlying population. A proper vector notation will give us the advantage of writing down linear transformations of a whole multidimensional vector distribution in a short and concise form.

Besides introducing random vectors and their components the next post

Multivariate Normal Distributions – II – random vectors and their covariance matrix

will also discuss related probability densities, expectation values and the definition of a covariance matrix for a random vector. Some simple properties of the covariance matrix will help us in further posts.

Affine properties of shear transformations

Properties of shear transformations and their matrices Msh

2D: Shearing of a square/rectangle

Shear transformation of a cube

Shearing of a cube: Surface layers

Remarks on a significant deficit of Matplotlib for 3D

Combinations of shear transformations

Combination of upper and lower triangular shear matrices?

By combining upper and lower triangular shear matrices we leave the group of pure shear operations

Just for fun: “Symmetric shearing” applied to a cube

Conclusion

Shear operations break symmetries of the affected figures

A shear operation preserves the orientation of tangent planes at some surface points

Coupling of more than two vector components – and a restriction

A first example produced with the help of Blender

Shear transformations are mathematically done by unipotent upper triangular matrices

Not a full and neither a symmetric matrix

Why did we not use a lower triangular matrix?

What about the diagonal elements of a shear matrix?

Basic properties of a shear matrix

A shear matrix cannot be diagonalized

More illustrations from Blender

Conclusion

Random Vectors and multivariate distributions

Samples

Random vectors

Linear operations on random vectors

Invertible transformations?

Probability density functions in the ℝn

Correlations of random vector components

Probability density after linear transformations and contour hyper-surfaces

Expectation value and covariance of a statistical vector distribution

Useful properties of the covariance matrix of a random vector

Remark on independence and correlation

Conclusion

Multivariate normal distribution and their projection into a 3-dimensional sub-space

Objective of this post

The covariance matrix determines everything

Creation of a MND from a standardized normal distribution

Creation of a trivariate normal distribution

Create some nested transparent ellipsoidal surfaces

SVD decomposition

Code snippets

Example views

Conclusion

ML is math in the end – and MNDs are no exception

Objectives of this post series

Objectives of the next post

Properties of shear transformations and their matrices M_sh

Probability density functions in the ℝⁿ