text:
# The Question

How does one prove that Rindler's definition of the covariant derivative of a covariant vector field $\lambda_a$ as \begin{align} \lambda_{a;c} = \lambda_{a,c} - \Gamma^{b}_{\ \ ca} \lambda_{b} \tag{Rindler} \end{align} cannot be written in terms of an eigendecomposition of the partial derivative such as \begin{align} \lambda_{a;c} = \dfrac{\partial y_d}{\partial x_a} \dfrac{\partial }{\partial x_c} \left( \dfrac{\partial x_b}{\partial y_d} \lambda_{b} \right) \tag{linuxfreebird} ? \end{align}

# The Definitions

The following symbols are defined by Rindler and interpreted by linuxfreebird:

\begin{align} \Gamma^{b}_{\ \ ca} = g^{bk} \Gamma_{k ca} \tag{Christoffel symbol} \end{align}

are the Christoffel symbols and

\begin{align} \tag{inverse metric tensor} g^{bd} = \dfrac{\partial x_b}{\partial y_k} \dfrac{\partial x_d}{\partial y_k} \end{align}

is the inverse of the metric tensor $g_{bd}$.

The eigendecomposition of a partial derivative is analogous to eigendecomposition of matrix operators. Given a matrix operator $A$, one can express $A$ in terms of its eigenvalue diagonal matrix $\Lambda$ and its eigenvector matrix $V$ such that $A = V\Lambda V^{-1}$.

# The Reason for the Question

Taken from Rindler, the Riemann tensor $R^{d}_{\ \ abc}$ is defined using the commutator of the covariant derivative from Eq. (Rindler) as

\begin{align} \lambda_{a;bc} - \lambda_{a;cb} = R^{d}_{\ \ abc} \lambda_{d} \tag{Riemman tensor}. \end{align}

However, if one attempts to calculate the Riemman tensor from Eq. (Riemman tensor) by using the definition of the covariant derivative from Eq. (linuxfreebird), one obtains the following result

\begin{align} \lambda_{a;bc} &=& \dfrac{\partial y_n}{\partial x_a} \dfrac{\partial }{\partial x_b} \left( \dfrac{\partial x_m}{\partial y_n} \dfrac{\partial y_f}{\partial x_m} \dfrac{\partial }{\partial x_c} \left( \dfrac{\partial x_d}{\partial y_f} \lambda_{d} \right) \right) \\&=& \dfrac{\partial y_n}{\partial x_a} \dfrac{\partial^2 }{\partial x_b \partial x_c} \left( \dfrac{\partial x_d}{\partial y_n} \lambda_{d} \right) , \end{align}

which states that the commutator of the covariant derivative is zero, hence the Riemman tensor is zero. This was the inconsistency I discovered from Eq. (linuxfreebird) and proves Eq. (linuxfreebird) cannot be an equivalent representation of Eq. (Rindler). However, I have been unable to prove why Eq. (linuxfreebird) is incorrect. Section (The Supporting Material) provides my formulation of Eq. (linuxfreebird).

# The Supporting Material

Below is my original formulation of Eq. (linuxfreebird):

Rindler defines the Christoffel symbol of the first kind as \begin{align} \tag{12} \Gamma_{b ca} = \dfrac{1}{2} \left( \dfrac{\partial g_{bc}}{\partial x_{a}} + \dfrac{\partial g_{ab}}{\partial x_{c}} - \dfrac{\partial g_{ca}}{\partial x_{b}} \right) , \end{align} where \begin{align} \tag{13} g_{ab} = \dfrac{\partial y_k}{\partial x_{a}} \dfrac{\partial y_k}{\partial x_{b}} \end{align} is the metric tensor. The partial derivative of $g_{ab}$ is \begin{align} \tag{14} \dfrac{\partial g_{ab}}{\partial x_{c}} = A_{acb} + A_{bca} , \end{align} where \begin{align} \tag{15} A_{acb} = \dfrac{\partial y_k}{\partial x_{a}} \dfrac{\partial^2 y_k}{\partial x_{c}\partial x_{b}} . \end{align} One can express $\Gamma_{bca}$ in terms of $A_{abc}$ as \begin{align} \tag{16} \Gamma_{b ca} = \dfrac{1}{2} \left( A_{bac} + A_{cab} + A_{acb} + A_{bca} - A_{cba} - A_{abc} \right) = A_{bac} , \end{align} which proves that \begin{align} \tag{17} \Gamma_{b ca} = \dfrac{\partial y_k}{\partial x_{b}} \dfrac{\partial^2 y_k}{\partial x_{c}\partial x_{a}} . \end{align} Rindler defines the Christoffel symbol of the second kind as \begin{align} \tag{18} \Gamma^{b}_{\ ca} = g^{bk} \Gamma_{k ca} , \end{align} where \begin{align} \tag{19} g^{ab} = \dfrac{\partial x_{a}}{\partial y_k} \dfrac{\partial x_{b}}{\partial y_k} \end{align} is the inverse of $g_{ab}$. One can express $\Gamma^{b}_{\ ca}$ in terms of the Jacobian components $\partial x_{a}/\partial y_{k}$ and prove the following equivalent expressions \begin{align} \tag{20} \Gamma^{b}_{\ ca} &=& \left( \dfrac{\partial x_{b}}{\partial y_n} \dfrac{\partial x_{k}}{\partial y_n} \right) \left( \dfrac{\partial y_m}{\partial x_{k}} \dfrac{\partial^2 y_m}{\partial x_{c}\partial x_{a}} \right) \\ &=& \tag{21} \dfrac{\partial x_{b}}{\partial y_n} \left( \dfrac{\partial x_{k}}{\partial y_n} \dfrac{\partial y_m}{\partial x_{k}} \right) \dfrac{\partial^2 y_m}{\partial x_{c}\partial x_{a}} \\ &=& \tag{22} \dfrac{\partial x_{b}}{\partial y_n} \left( \delta_{mn} \right) \dfrac{\partial^2 y_m}{\partial x_{c}\partial x_{a}} \\ &=& \tag{23} \dfrac{\partial x_{b}}{\partial y_m} \dfrac{\partial }{\partial x_{c}} \left( \dfrac{\partial y_m}{\partial x_{a}} \right) \\ &=& \tag{24} \dfrac{\partial }{\partial x_{c}} \left( \dfrac{\partial x_{b}}{\partial y_m} \dfrac{\partial y_m}{\partial x_{a}} \right) - \dfrac{\partial y_m}{\partial x_{a}} \dfrac{\partial }{\partial x_{c}} \left( \dfrac{\partial x_{b}}{\partial y_m} \right) \\ &=& \tag{25} \dfrac{\partial }{\partial x_{c}} \left( \delta_{ab} \right) - \dfrac{\partial y_m}{\partial x_{a}} \dfrac{\partial }{\partial x_{c}} \left( \dfrac{\partial x_{b}}{\partial y_m} \right) \\ &=& \tag{26} - \dfrac{\partial y_m}{\partial x_{a}} \dfrac{\partial }{\partial x_{c}} \left( \dfrac{\partial x_{b}}{\partial y_m} \right) . \end{align} Rindler defines the covariant derivative of a covariant vector field $\lambda_{a}$ as \begin{align} \tag{27} \lambda_{a;c} = \lambda_{a,c} - \Gamma^{b}_{\ ca} \lambda_{b} . \end{align} One can define the covariant derivative in terms of $\partial x_{a}/\partial y_{k}$ as \begin{align} \tag{28} \lambda_{a;c} &=& \dfrac{\partial \lambda_{a} }{\partial x_{c}} + \dfrac{\partial y_m}{\partial x_{a}} \dfrac{\partial }{\partial x_{c}} \left( \dfrac{\partial x_{b}}{\partial y_m} \right) \lambda_{b} \\&=& \tag{29} (\delta_{ab}) \dfrac{\partial \lambda_{b} }{\partial x_{c}} + \dfrac{\partial y_m}{\partial x_{a}} \dfrac{\partial }{\partial x_{c}} \left( \dfrac{\partial x_{b}}{\partial y_m} \right) \lambda_{b} \\&=& \tag{30} \left( \dfrac{\partial y_m}{\partial x_{a}} \dfrac{\partial x_{b}}{\partial y_m} \right) \dfrac{\partial \lambda_{b} }{\partial x_{c}} + \dfrac{\partial y_m}{\partial x_{a}} \dfrac{\partial }{\partial x_{c}} \left( \dfrac{\partial x_{b}}{\partial y_m} \right) \lambda_{b} \\&=& \tag{31} \dfrac{\partial y_m}{\partial x_{a}} \left( \dfrac{\partial x_{b}}{\partial y_m} \dfrac{\partial (\lambda_{b} ) }{\partial x_{c}} + \dfrac{\partial }{\partial x_{c}} \left( \dfrac{\partial x_{b}}{\partial y_m} \right) \lambda_{b} \right) \\&=& \tag{32} \dfrac{\partial y_m}{\partial x_{a}} \dfrac{\partial }{\partial x_{c}} \left( \dfrac{\partial x_{b}}{\partial y_m} \lambda_{b} \right) . \end{align}

-
The indices in your equation 5 make no sense... – Danu Mar 1 '14 at 19:34
Also, could it be that your equation 4 is just saying that $\nabla_a\lambda_b=\frac{1}{\sqrt{-|g|}}\partial_a(\sqrt{-|g|}\lambda_b)$? Because that much is true – Danu Mar 1 '14 at 19:36
It would help if you provided a specific reference to where you found the alternative definition. – Emilio Pisanty Mar 2 '14 at 21:54
A bunch of your equations here are not inverting the index of the denominator ($\frac1{B_\mu} = C^{\mu}$, not $C_\mu$ for some quantities $B,C$) – Manishearth Mar 3 '14 at 3:38
While more information never hurts, there is still quite a bit you could do to make your question more readable and understandable. It should be clear without scrolling down what the question is about, and where the 'supplementary details' begin. Don't be afraid to use extra formatting like **bold** or separators like ---. – Emilio Pisanty Mar 4 '14 at 9:04

Let there be given a manifold $(M,\nabla)$ equipped with a (not necessarily torsionfree) tangent bundle connection $\nabla$.

I got the (possibly faulty) impression from reading the first lines in OP's question formulation (v18) that OP is asking:

Is it possible that the local coordinate expression for the covariant derivative of a co-vector/one-form $\lambda =\lambda_a \mathrm{d}x^a$ (in some local coordinate system $x^a$) could be on the form $$\tag{1} (\nabla_c \lambda)_a~=~\frac{\partial y^d}{\partial x^a} \frac{\partial }{\partial x^c}\left(\frac{\partial x^b}{\partial y^d} \lambda_b \right)$$ for some (invertible, smooth, locally defined) functions $y^d=y^d(x)$?

Is it possible that the Christoffel symbols$^1$ (in some local coordinate system $x^a$) could be of the form $$\tag{2} \Gamma^{(x)b}_{ca}~=~\frac{\partial^2 y^d}{\partial x^c\partial x^a} \frac{\partial x^b}{\partial y^d}$$ for some (invertible smooth, smooth, locally defined) functions $y^d=y^d(x)$?

Recall that the Christoffel symbol $\Gamma^{\lambda}_{\mu\nu}$ does not transform as a tensor under a local coordinate transformation $x^a \to y^b=y^b(x)$, but rather with an inhomogeneous term, which is built from the second derivative of the coordinate transformation

$$\tag{3}\frac{\partial y^f}{\partial x^c} \Gamma^{(x)c}_{ab} ~=~\frac{\partial y^d}{\partial x^a}\, \frac{\partial y^e}{\partial x^b}\, \Gamma^{(y)f}_{de} +\frac{\partial^2 y^f}{\partial x^a \partial x^b}.$$

Therefore formulation (2) is equivalent of asking:

Does there exists a local coordinate system $y^d$ such the Christoffel symbols vanish $$\tag{4} \Gamma^{(y)c}_{ab}~=~0?$$

This, in turn, is equivalent to asking

Is $(M,\nabla)$ flat?

Does there exist a local coordinate system $y^d$ such that the metric tensor $$\tag{5} g^{(y)}_{ab}~=~\eta_{ab}$$ is constant?

In an arbitrary coordinate system $x^a$, is the metric tensor of the form $$\tag{6} g^{(x)}_{ab}~=~\frac{\partial y^c}{\partial x^a}\eta_{cd} \frac{\partial y^d}{\partial x^b}$$ for some (invertible smooth, smooth, locally defined) functions $y^d=y^d(x)$?

Rereading OP's question formulation, condition (6) seems to be one of OP's assumptions. If we are granted assumption (6), then OP's initial claim (1) holds.

--

$^1$ It is convenient to call $\Gamma^c_{ab}$ Christoffel symbols even if the tangent-space connection $\nabla$ is not torsionfree.

-
@linuxfreebird: I updated the answer. Note that eq. numbers may have shifted. – Qmechanic Mar 6 '14 at 14:33

Equation (13) expresses the metric on an embedded hypersurface given by the relations $y^k = y^k(x^a)$. However, the equation for the inverse metric (4-th equation) is in general not correct.

Take for example a hypersurface defined by: $y^1 = x^1$, $y^2 = x^2$, $y^3 = x^2$. In our case, the partial derivative of $x^2$ with respect to $y^2$ or $y^3$ cannot be defined. For example, in order to differentiate with respect to $y^2$ we must vary $y^2$ while keeping $y^3$ constant. This is impossible to do on the hypersurface.

Update: This is a detailed answer to linuxfreebird comment.

Although I did not receive my education from Rindler's book, the material of the subject of this question appears in many pedagogical introductions of general relativity, in which an $N$ dimensional space time $X$ is isometrically embedded in a larger flat $M$-dimensional flat space $Y$:

$$\begin{matrix} X \rightarrow Y \\ y^k = y^k(x^1, ...., x^N), k = 1, ..., M \end{matrix}$$

The embedding is given as a set of equations defining $X$ as a hypersurface into $Y$. The functions need not be invertible as we want them to describe embeddings of spaces of different dimensions.

Equation (4) is just the relation expressing that the embedding is isometric. This can be seen more explicitly if we write this equation in the form:

$$g_{a b} =\eta_{kl} \frac{\partial y^k }{\partial x^a} \frac{\partial y^l}{\partial x^b}$$

Where $\eta$ is the flat metric on the manifold $Y$. Now, this equation looks as a metric transformation, i.e., it preserves distances.

This approach is called extrinsic geometry. It can be used to define the various geometrical objects of the space-time $X$ (such as the Levi-Civita connection and the Riemann curvature) in terms of the embedding equations. This approach is correct, because for any space-time $X$ there exists such an embedding for some high enough $M$ (This is called Nash embedding theorem), but it is seldom used in general relativity real work. Most work in general relativity uses intrinsic geometry without resorting to any embedding equations.

The key point is to remember that the $y$s are functions of the $x$s and every time a partial derivative of $x$ with respect to $y$ is used, means that these equations were inverted, but this cannot be done in general, except in the particular case where these equations are invertible.

Now, in the case when the equations are invertible (it necessarily implies that $M=N$ ), everything is correct, but this case refers only to space-times that can be isometrically embedded into a flat space-time of the same dimension. These space-times are necessarily flat, thus not general enough to cover interesting cases in general relativity.

May be you can profit from this analysis if you continue your exercise and prove that the Riemann tensor computed from your form of the Levi-Civita connection is globally identically zero.

-
@linuxfreebird. In the case when the transformation is invertible, then your analysis is correct, but this special case is not really interesting in general relativity. I updated my answer with a detailed explanation. – David Bar Moshe Mar 6 '14 at 13:15
Nash's theorem is for the Riemannian case. Does apply for Lorenzian manifold as well? My guess is that there is an analogous result, do you know any references? – MBN Mar 6 '14 at 13:35
@MBN there are partial results, for example Ellis meghnad.iucaa.ernet.in/~tarun/pprnt/topology/… mentions a result by Clarke that every 3+1 noncompact space-time can be isometrically embedded into a flat $\mathbb{R}^{87,3}$. – David Bar Moshe Mar 6 '14 at 14:14
@DavidBarMoshe: Thanks, that looks interesting. – MBN Mar 6 '14 at 14:42

Special thanks to Qmechanic and David Bar Moshe for providing their answers for this post. Combining Qmechanic's and David Bar Moshe's insights, I finally figured out the confusing element in the problem.

The covariant derivative can be written in terms of an eigendecomposition of the partial derivative, the Christoffel symbols can be written in terms of non-invertible extensive variables, and the Riemman tensor is not necessarily equal to zero. Here is my answer:

I have written the covariant derivative \begin{align} \lambda_{a;c} = \dfrac{\partial y^k}{\partial x^a} \dfrac{\partial }{\partial x^c} \left( \dfrac{\partial y_k}{\partial x^p} g^{pb} \lambda_{b} \right) \end{align} in terms of its eigendecomposition and indices "a" and "c" and again \begin{align} \lambda_{b;d} = \dfrac{\partial y^j}{\partial x^b} \dfrac{\partial }{\partial x^d} \left( \dfrac{\partial y_j}{\partial x^q} g^{qn} \lambda_{n} \right) \end{align} with different indices "b" and "d". Combining the covariant derivatives forms the following: \begin{align} \lambda_{a;cd} =\\ \dfrac{\partial y^k}{\partial x^a} \dfrac{\partial }{\partial x^c} \left( \Lambda^{\ \ j}_{k} \dfrac{\partial }{\partial x^d} \left( \dfrac{\partial y_j}{\partial x^q} g^{qn} \lambda_{n} \right) \right) \end{align} where \begin{align} \Delta^{\ \ j}_{k} = \dfrac{\partial y_k}{\partial x^p} g^{pb} \dfrac{\partial y^j}{\partial x^b} \neq \delta_{k}^{\ \ j} \end{align} because \begin{align} g^{pb} \neq \dfrac{\partial x^p}{\partial y^n} \dfrac{\partial x^b}{\partial y_n} \end{align}

is ill-defined.

The Riemann tensor $R_{acd}^{ \ \ \ \ m}$ is created by computing the commutator of the covariant derviatives as follows:

\begin{align} \lambda_{a;[c,d]} = \left( \dfrac{\partial y^k}{\partial x^a} \dfrac{\partial \Delta^{\ \ j}_{k} }{\partial x^c} \dfrac{\partial }{\partial x^d} \left( \dfrac{\partial y_j}{\partial x^p} g^{pm} \right) \right)_{[c,d]} = R_{acd}^{ \ \ \ \ m} \lambda_m . \end{align}

If the number of dimensions of $y$ equals the number of dimensions in $x$, then $\Delta^{\ \ j}_{k} = \delta^{\ \ j}_{k}$, thus $R_{acd}^{ \ \ \ \ m} = 0$.

-


text:
# Rotate 3d plane

I have a plane in 3D space that formed from 3 poin $P_1=(x_1, y_1, z_1)$, $P_2=(x_2, y_2, z_2)$, $P_3=(x_3, y_3, z_3)$

I want to rotate and transform this points (equally related plane) into 2D space (Avoiding $z$ axis but save distance and relations in 2D plane).

• To clarify, you want to rotate the plane containing the points $P_1,\,P_2,\,P_3$ into the $xy$ plane. – Daryl Apr 29 '13 at 8:08
• It's only necessary to rotate the plane into a plane parallel to the x-y plane. (If the plane has any distance from the origin, then it will be impossible to rotate it into the x=y plane, which passes thru the origin). Once rotated parallel to x-y plane, all the z-coordinates in the original plane will be the same, and can be ignored, making it a 2D plane. – pbierre Aug 2 '17 at 2:59

As a first step, I'd move $P_1$ to the origin, so that the points become $P_1=(0,0,0)$, $P_2=(x_2-x_1, y_2-y_1, z_2-z_1)$, $P_3=(x_3-x_1, y_3-y_1, z_3-z_1)$. From there it's just a question of applying a rotation matrix.

A rotation matrix $\mathbf{R}$ which preserves all distances and shapes etc is given by

$$\mathbf{R}=\left( \begin{array} [c]{ccc}% e_{1}^{x} & e_{2}^{x} & e_{3}^{x}\\ e_{1}^{y} & e_{2}^{y} & e_{3}^{y}\\ e_{1}^{z} & e_{2}^{z} & e_{3}^{z}% \end{array} \right)$$

where the vectors $\left( e_{1}^{\alpha},e_{2}^{\alpha},e_{3}^{\alpha }\right)$ for $\alpha\in\left\{ x,y,z\right\}$ are orthogonal and of unit length, i.e.

$$\sum_{i=1}^{3}e_{i}^{\alpha}e_{i}^{\beta}=\left\{ \begin{array} [c]{ccc}% 1 & & \alpha=\beta\\ 0 & & \alpha\neq\beta \end{array} \right. .%$$

Applying $\mathbf{R}$, your rotated points $P_{i}^{\prime}=\left( x_{i}^{\prime},y_{i}^{\prime},z_{i}^{\prime}\right)$ are given by

$$\left( \begin{array} [c]{c}% x_{i}^{\prime}\\ y_{i}^{\prime}\\ z_{i}^{\prime}% \end{array} \right) =\mathbf{R}\left( \begin{array} [c]{c}% x_{i}\\ y_{i}\\ z_{i}% \end{array} \right)$$

I suggest making the top row of $\mathbf{R}$ the unit vector in the direction from $P_{1}$ to $P_{2}$, i.e.

$$\left( \begin{array} [c]{c}% e_{1}^{x}\\ e_{2}^{x}\\ e_{3}^{x}% \end{array} \right) =\frac{1}{\sqrt{\left( x_{1}-x_{2}\right) ^{2}+\left( y_{1}% -y_{2}\right) ^{2}+\left( z_{1}-z_{2}\right) ^{2}}}\left( \begin{array} [c]{c}% x_{1}-x_{2}\\ y_{1}-y_{2}\\ z_{1}-z_{2}% \end{array} \right)$$

which would mean that the line from $P_{1}$ to $P_{2}$ gets mapped into the $x$-axis.

Working out what you could use for the second and third rows of $\mathbf{R}$ is now just a question of solving some simple linear equations, to make sure that the line from $P_{1}$ to $P_{3}$ has no z component.

To solve for $\left( e_{1}^{z},e_{2}^{z},e_{3}^{z}\right)$, the vector $e_{i}^{z}$ must have a zero dot-product with both $e_{i}^{x}$ and $\left( x_{3}-x_{1},y_{3}-y_{1},z_{3}-z_{1}\right)$, so the equations are

\begin{align*} e_{1}^{x}e_{1}^{z}+e_{2}^{x}e_{2}^{z}+e_{3}^{x}e_{3}^{z} & =0\\ \left( x_{3}-x_{1}\right) e_{1}^{z}+\left( y_{3}-y_{1}\right) e_{2}% ^{z}+\left( z_{3}-z_{1}\right) e_{3}^{z} & =0 \end{align*}

Hence

$$\left( \begin{array} [c]{c}% e_{1}^{z}\\ e_{2}^{z}\\ e_{3}^{z}% \end{array} \right) =\lambda_{z}\left( \begin{array} [c]{c}% \left( z_{3}-z_{1}\right) e_{2}^{x}-\left( y_{3}-y_{1}\right) e_{3}^{x}\\ \left( x_{3}-x_{1}\right) e_{3}^{x}-\left( z_{3}-z_{1}\right) e_{1}^{x}\\ \left( y_{3}-y_{1}\right) e_{1}^{x}-\left( x_{3}-x_{1}\right) e_{2}^{x}% \end{array} \right)$$

for some $\lambda_z$, which should be determined so that $e_{i}^{z}$ is a vector of unit length. Finally $e_{i}^{y}$ can be determined as the vector which is perpendicular to both $e_{i}^{x}$ and $e_{i}^{z}$, i.e.

$$\left( \begin{array} [c]{c}% e_{1}^{y}\\ e_{2}^{y}\\ e_{3}^{y}% \end{array} \right) =\lambda_{y}\left( \begin{array} [c]{c}% e_{3}^{z}e_{2}^{x}-e_{2}^{z}e_{3}^{x}\\ e_{1}^{z}e_{3}^{x}-e_{3}^{z}e_{1}^{x}\\ e_{2}^{z}e_{1}^{x}-e_{1}^{z}e_{2}^{x}% \end{array} \right)$$

where again $\lambda_{y}$ is determined so that $e_{i}^{y}$ is a vector of unit length.

As the accepted answer suggests, first apply a translation of $-P_1$ to all points. Your three points are then $O,$ $A,$ and $B$ (where $A=P_2-P_1$ and $B=P_3-P_1$). Now compute the plane unit normal: $$n = \frac{A \times B}{| A\times B|}.$$

To rotate the plane to the $xy$ plane, all you need to do is rotate $n$ to $(0,0,1)$. The simplest way to do this is to rotate $n$ about the $z$ axis into the positive $x$ half of the $xz$ plane and then rotate about the $y$ axis until $n$ is parallel to the positive $z$ axis.

You could also use Rodrigues' rotation formula to rotate about the axis $n\times (0,0,1),$ but this rotation matrix is much more complicated than the matrix for a rotation about a coordinate axis.

Computing the rotation matrix that rotates $n$ about the $z$ axis into the positive $x$ half of the $xz$ plane does not require computing an angle. The rotation matrix is simply given in terms of $n$: $$R_z = \left( \begin{array}{lll} n_x/\sqrt{n_x^2+n_y^2} & n_y/\sqrt{n_x^2+n_y^2} & 0\\ -n_y/\sqrt{n_x^2+n_y^2} & n_x/\sqrt{n_x^2+n_y^2} & 0\\ 0 & 0 & 1 \end{array} \right).$$ The rotation matrix that rotates $n'=R_zn$ about the $y$ axis to $(0,0,1)$ also does not require computing an angle. It is $$R_y = \left( \begin{array}{lll} n'_z & 0 & -n'_x \\ 0 & 1 & 0\\ n'_x & 0 & n'_z \end{array} \right).$$

You sound like you could benefit from learning how to represent direction in 3D using a direction vector. This entails learning how to normalize a vector down to unit length, which results in a direction vector.

Then, once comfy with that, learn how to solve for the orientation of (direction perpendicular to) the plane of a 3D triangle:

  o <-- normalize[(P2-P1) x (P3-p2)]   (vector cross product of 2 vector diffs)

Now, you want to define the desired 3D coordinate rotation which aligns the orientation of your triangle-plane with the z-axis. Again, using cross products, these are the 3 columns of the 3x3 rotation matrix you want:

R = [ newXaxis newYaxis newZaxis ] (three orthogonal direction vectors)

By choice, you choose your value of o as newZaxis.

The other two new axes have to be orthogonal to newZaxis and to each other. The normalized vector cross product gives you the mutual orthogonal direction to any two distinct directions d1 an d2 (so long as they are not coincident or opposite directions):

newXaxis = normalize(newZaxis x [ 1 0 0 ] )

(if newZaxis == [ 1 0 0 ], use newXaxis = normalize(newZaxis x [ 0 1 0 ] ) instead

newYaxis = newZaxis x newYaxis

Now, coordinate rotate each point of your triangle by rotator R:

p1' <-- R • p1 (vector times a matrix)

p2' <-- R • p2

p3' <-- R • p3

You'll notice that p1', p2' and p3' all have identical z-coords. You can throw these out, and just use their x and y coordinates. You've effectively projected your triangle into a 2D x-y plane.

format: plain


text:
• 12
• 14
• 13
• 10
• 11
• ### Similar Content

• By ucfchuck
I am feeding in 16 bit unsigned integer data to process in a compute shader and i need to get a standard deviation.
So I read in a series of samples and push them into float arrays
float vals1[9], vals2[9], vals3[9], vals4[9]; int x = 0,y=0; for ( x = 0; x < 3; x++) { for (y = 0; y < 3; y++) { vals1[3 * x + y] = (float) (asuint(Input1[threadID.xy + int2(x - 1, y - 1)].x)); vals2[3 * x + y] = (float) (asuint(Input2[threadID.xy + int2(x - 1, y - 1)].x)); vals3[3 * x + y] = (float) (asuint(Input3[threadID.xy + int2(x - 1, y - 1)].x)); vals4[3 * x + y] = (float) (asuint(Input4[threadID.xy + int2(x - 1, y - 1)].x)); } } I can send these values out directly and the data is as expected

Output1[threadID.xy] = (uint) (vals1[4] ); Output2[threadID.xy] = (uint) (vals2[4] ); Output3[threadID.xy] = (uint) (vals3[4] ); Output4[threadID.xy] = (uint) (vals4[4] ); however if i do anything to that data it is destroyed.
vals1[4] = vals1[4]/2;
or a
vals1[4] = vals[1]-vals[4];
the data is gone and everything comes back 0.

How does one go about converting a uint to a float and performing operations on it and then converting back to a rounded uint?
• By fs1
I have been trying to see how the ID3DInclude, and how its methods Open and Close work.
I would like to add a custom path for the D3DCompile function to search for some of my includes.
I have not found any working example. Could someone point me on how to implement these functions? I would like D3DCompile to look at a custom C:\Folder path for some of the include files.
Thanks
• By stale
I'm continuing to learn more about terrain rendering, and so far I've managed to load in a heightmap and render it as a tessellated wireframe (following Frank Luna's DX11 book). However, I'm getting some really weird behavior where a large section of the wireframe is being rendered with a yellow color, even though my pixel shader is hard coded to output white.

The parts of the mesh that are discolored changes as well, as pictured below (mesh is being clipped by far plane).

Here is my pixel shader. As mentioned, I simply hard code it to output white:
float PS(DOUT pin) : SV_Target { return float4(1.0f, 1.0f, 1.0f, 1.0f); } I'm completely lost on what could be causing this, so any help in the right direction would be greatly appreciated. If I can help by providing more information please let me know.

• Hello,
i try to implement voxel cone tracing in my game engine.
At first step i try to emplement the easiest "poor mans" method
a.  my test scene "Sponza Atrium" is voxelized completetly in a static voxel grid 128^3 ( structured buffer contains albedo)
b. i dont care about "conservative rasterization" and dont use any sparse voxel access structure
c. every voxel does have the same color for every side ( top, bottom, front .. )
d.  one directional light injects light to the voxels ( another stuctured buffer )
I will try to say what i think is correct ( please correct me )
GI lighting a given vertecie  in a ideal method
A.  we would shoot many ( e.g. 1000 ) rays in the half hemisphere which is oriented according to the normal of that vertecie
B.  we would take into account every occluder ( which is very much work load) and sample the color from the hit point.
C. according to the angle between ray and the vertecie normal we would weigth ( cosin ) the color and sum up all samples and devide by the count of rays
Voxel GI lighting
In priciple we want to do the same thing with our voxel structure.
Even if we would know where the correct hit points of the vertecie are we would have the task to calculate the weighted sum of many voxels.
Saving time for weighted summing up of colors of each voxel
To save the time for weighted summing up of colors of each voxel we build bricks or clusters.
Every 8 neigbour voxels make a "cluster voxel" of level 1, ( this is done recursively for many levels ).
The color of a side of a "cluster voxel" is the average of the colors of the four containing voxels sides with the same orientation.

After having done this we can sample the far away parts just by sampling the coresponding "cluster voxel with the coresponding level" and get the summed up color.
Actually this process is done be mip mapping a texture that contains the colors of the voxels which places the color of the neighbouring voxels also near by in the texture.
Cone tracing, howto ??
Here my understanding is confus ?? How is the voxel structure efficiently traced.
I simply cannot understand how the occlusion problem is fastly solved so that we know which single voxel or "cluster voxel" of which level we have to sample.
Supposed,  i am in a dark room that is filled with many boxes of different kind of sizes an i have a pocket lamp e.g. with a pyramid formed light cone
- i would see some single voxels near or far
- i would also see many different kind of boxes "clustered voxels" of different sizes which are partly occluded
How do i make a weighted sum of this ligting area ??
e.g. if i want to sample a "clustered voxel level 4" i have to take into account how much per cent of the area of this "clustered voxel" is occluded.
Please be patient with me, i really try to understand but maybe i need some more explanation than others
best regards evelyn

• Hi guys, when I do picking followed by ray-plane intersection the results are all wrong. I am pretty sure my ray-plane intersection is correct so I'll just show the picking part. Please take a look:

// get projection_matrix DirectX::XMFLOAT4X4 mat; DirectX::XMStoreFloat4x4(&mat, projection_matrix); float2 v; v.x = (((2.0f * (float)mouse_x) / (float)screen_width) - 1.0f) / mat._11; v.y = -(((2.0f * (float)mouse_y) / (float)screen_height) - 1.0f) / mat._22; // get inverse of view_matrix DirectX::XMMATRIX inv_view = DirectX::XMMatrixInverse(nullptr, view_matrix); DirectX::XMStoreFloat4x4(&mat, inv_view); // create ray origin (camera position) float3 ray_origin; ray_origin.x = mat._41; ray_origin.y = mat._42; ray_origin.z = mat._43; // create ray direction float3 ray_dir; ray_dir.x = v.x * mat._11 + v.y * mat._21 + mat._31; ray_dir.y = v.x * mat._12 + v.y * mat._22 + mat._32; ray_dir.z = v.x * mat._13 + v.y * mat._23 + mat._33;
That should give me a ray origin and direction in world space but when I do the ray-plane intersection the results are all wrong.
If I click on the bottom half of the screen ray_dir.z becomes negative (more so as I click lower). I don't understand how that can be, shouldn't it always be pointing down the z-axis ?
I had this working in the past but I can't find my old code

# DX11 Directxtex; Unresolved External, Lib Needed?

This topic is 621 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

Hi all,

So far I thought I could just download the .H and .CPP files, place them in a folder which I include in my project, and go for it.

Unfortunaly that doesn't work, my solution/ code compiles, but I get an unresolved external (like a LIB is missing or something).

Do you have any experience with using this library (in particular the DDS loader)?
On github I also see there's a Visual Studio project of the DDS loader/ directxtex library, I could compile that and generate a LIB, but I don't know if that's really needed to just use the DDS loader.

Any input is appreciated.

#include <DDSTextureLoader.h>

// load and create 1 texture
//HR(D3DX11CreateShaderResourceViewFromFile(md3dDevice, L"Textures/darkbrickdxt1.dds", 0, 0, &mDiffuseMapSRV, 0 ));

HR(CreateDDSTextureFromFile(md3dDevice, L"Textures/darkbrickdxt1.dds", NULL, &mDiffuseMapSRV));

1>main.obj : error LNK2001: unresolved external symbol "long __cdecl DirectX::CreateDDSTextureFromFile(struct ID3D11Device *,wchar_t const *,struct ID3D11Resource * *,struct ID3D11ShaderResourceView * *,unsigned int,enum DirectX::DDS_ALPHA_MODE *)" (?CreateDDSTextureFromFile@DirectX@@YAJPAUID3D11Device@@PB_WPAPAUID3D11Resource@@PAPAUID3D11ShaderResourceView@@IPAW4DDS_ALPHA_MODE@1@@Z)
1>E:\projects\D3D11 practice\Minimal application\Release\Minimal application.exe : fatal error LNK1120: 1 unresolved externals

Edited by cozzie

##### Share on other sites

error LNK2001: unresolved external symbol while looking for CreateDDSTextureFromFile. Did you forget to add 'DDSTextureLoader.cpp' to your project? :)

A build of DirectXTex.lib from DirectXTex-master is not required for the light-weight DDS runtime loader.

##### Share on other sites
Thanks. I was afraid I had to do that. Was hoping that I could use it as an "external dependency" instead of add the source to my projects.

On github I saw that the dds and wic loader are not part of the directxtex project, so compiling that to a lib wont work either.

I'll just have to manage the files within my project.

##### Share on other sites

You could just create a separate project as static lib for those files yourself.


text:
# There are enough points, but how do you determine the equation?

Algebra Level 5

There is a hyberbola that passes through points

$A = (0, 6)$

$B = (1, 3)$

$C = (3, 1)$

$D = (6, 0)$

$E = (10, -1)$

If the equation for this hyperbola can be represented as

$ax^2 + by^2 + cxy + dx + ey + f = 0$

where the $$\gcd (a, b, c, d, e, f)$$ is 1, find the minimum value of

$a + b + c + d + e + f$

×


text:
19 Pages

## Introduction

Statistics are important communication tools! They're a way for you to tell the reader how confident you are in your measurements and conclusions, and a way for you to tell us that you know how dodgy reality can be. Here, two topics are of importance: Propagation of Error and Regression. Since you've all taken your stats course already some details will be omitted.

A simple pump calibration experiment will be described to give us a data set to work with, and MATLAB will be used for all analyses and plotting.

## Reporting the Error

### Sources of Error

If you report an experimentally-determined value then it should have an uncertainty (error) to go along with it. We typically report error in one of two ways: as a number or a plotted point. The former requires error to be reported explicitly (as in 23 +/- 2 mL/s) whereas the latter uses error bars to graphically indicate the error.

Whenever you report an error you should also include a brief description of how the error was determined, but you need do this only once per error term. For example, if you're reporting tank concentration measurements and you had to use a calibration curve then the first measurement you report should be something like, "The tank concentration was 100 +/- 8 ppm, where the uncertainty is the prediction interval of the equipment calibration curve." Subsequent reporting of concentrations made in the same way can be reported without the explanation.

There are four common sources of error in our lab experiments:

1. Measured variables. If you use a piece of equipment to directly measure a variable then you can usually find an error estimate on the equipment itself. For example, graduated cylinders usually have the words "Volume (mL) +/- 5%" printed right on the cylinder. If you're using equipment that doesn't provide uncertainty then the error is half of the smallest quantified value (e.g., if a rule has marks every 1 mm then the error is 0.5 mm, or if a scale reads to three digits of precision (0.000) then the error is 0.0005). If a quantity is calculated from several measured variables then you'll need to perform propagation of error (see below) to report the total error. In both cases, the first such error reported would include a descriptive phrase such as, "The flow rate was 5 +/- 1 mL/s, where the uncertainty results from propagation of equipment uncertainties."
2. Coefficients from calibration curves. If you use a piece of equipment to directly measure a variable but you have to convert the output signal of the equipment to a meaningful quantity by means of a calibration curve then you should report the error as the prediction interval of the calibration curve. For example, if you measured the salinity of a tank by means of a conductivity probe then the probe output will have been in uS/cm (or something similar) and you'll have needed to convert to a meaningful concentration such as ppm by means of a calibration curve. In this case you'd report a measured salt concentration as something like, "The tank salt concentration was 1050 +/- 20 ppm where the uncertainty is the prediction interval from the equipment calibration curve."
3. Repeated measurements. If you make three or more independent measurements of a variable then you can calculate a 95% confidence interval[1] on the mean and use that as your error estimate. Confidence intervals are preferred over standard deviations because the confidence interval indicates the probability that the reported interval contains the true value of the parameter (usually a mean or regression coefficient), whereas the standard deviation only describes the distribution of the sampled data (which may or may not be a good estimate of the distribution of the true data). For example, if you repeat the synthesis of liposome nanoparticles three times under the same conditions (and therefore expect the same particle size) then you could say something like, "The mean particle diameter was 160 +/- 30 nm where the uncertainty is the 95% confidence interval of the mean."
4. Coefficients from regression analysis. In most experiments you'll produce sets of data for which theoretical correlations exist, such as the diffusion constants in the RO experiment or the rate constants in the UV Photocatalysis experiment. In these cases the coefficients that result from a regression analysis of the data are the parameters of interest and the error should be reported as the 95% confidence interval of the regression coefficient. In cases where the slope is related but not directly equal to the parameter of interest (as is the case with some RO or UV coefficients) then you'll have to apply propagation of error rules as well to "extract" the desired value. For example, you might report one of the empirical coefficients determined in the Fuel Cell experiment (see Eq. (6)) in the following manner: "The open-circuit cell potential parameter was determined to be $$E_{0,R} = 0.88 \pm 0.6 \textrm{ V}$$, where the uncertainty is the 95% confidence interval on the regression coefficient."

The material below describes methods to determine each of these error estimates, with the exception of confidence interval on the mean from repeated measurements which is discussed here (and is usually included in the output of any descriptive statistics function).

### Numerical Precision

For the purposes of CENG 176A/B reports, error should always be reported as $$\textrm{(measured value)}\pm\textrm{(error)}$$, where the error has the same units as the measure value. For simplicity, we'll use the following (conservative) rules regarding significant figures:

1. Round the error up to one significant figure.
2. Round and report the measured variable to the same precision (i.e., number of decimal points) as the error.

For example, if the measured value was 22.6654 mL/s and the uncertainty was 1.66 mL/s, we'd first round the uncertainty up to 2 mL/s (one significant figure) and then round the measured value to 23 mL/s (same number of decimal points). The measure variable would therefore be reported as 23 +/- 2 mL/s.

This is a good general procedure when procedures like regression or propagation of error are involved but it has a few limitations to be aware of (as do most rules for assigning significant figures):

1. If the error was produced from an instrument for which you're going to estimate error as half the smallest measurement (like a ruler or balance) then you report the measured value to the precision of the instrument along with half the smallest measurement as the uncertainty (in other words, don't round anything). For example, if you're using a balance with four digits of precision (e.g., 0.0000 g) and you measure a sample of 0.0221 g, then you should report this number as 0.0221 +/- 0.00005 g even though that's four digits of precision on the measured value and five digits of precision on the error (normally they'll match). Depending on which reference you read you'll see variations of this kind of rule, which can be annoying.
2. Occasionally this rule can produce very large errors but that's less of a problem and more of a feature. For example, suppose you determine by regression a coefficient of $$b = 2.4359 \pm 0.6331$$: according to the rounding rules the error should round to one significant figure ($$1$$) and the value to equal precision $$2$$ which results in $$b = 2.4 \pm 0.7$$. There's not really a good way around this (at least, none that I'm aware of) to determine in an unambiguous manner how many digits to retain (other than to use your engineering judgment to decide if more digits are warranted).
3. This method doesn't work for cases like the LPCVD simulation and liposome nanoparticle analysis because we can't see inside the "black box" of how the machine (or computer program) is calculating its numbers. As with the previous case we simply use our best engineering judgment to report a reasonable number of significant figures.

## The Experiment

Calibration curves are required whenever two variables are related but in a way that's too complex or difficult to theoretically derive. For example, consider an inexpensive, variable-speed pump that has a knob to adjust the flow rate, but the knob has marks of 1, 2, 3...5 with no units. To translate from these arbitrary (also called "nominal") values to a meaningful flow unit, we perform the following experiment:

1. Connect the pump to a reservoir so that water is pumped in a loop.
2. Set the pump knob to 1 and let the system stabilize.
3. Use a stopwatch and graduated cylinder to measure the volume discharged in 10 s.
4. Repeat 2-3 for all knob settings.

After completing Step 4, we assemble the data as shown in Table 1. We'd now like to derive an equation to allow us to easily convert from nominal pump markings to volumetric flow rate. Towards that end we'll use statistical analyses to determine both the equation and our confidence in the equation.

Table 1. Experimental data for the hypothetical pump calibration experiment.
Knob Setting Measured time (s) Measured volume (mL)
0 10.1 0
1 10.0 23
2 10.2 43
3 9.9 71
4 9.9 83
5 10.1 120

## Propagation of Error

The general formula for the error $$e_Y$$ of a quantity $$Y$$ that cannot be measured directly but is a function of one or more directly measured quantities $$y_i$$, each with their own error $$e_i$$, is

$e_Y^2 = \sum_{i=1}^{N} \left( \frac{\partial Y}{\partial y_i}\right)^2 e_i^2,$

where $$N$$ is the number of measured quantities, and the errors $$e_i$$ were assumed to be small and independent.[2]

### Example 1: Derivation of Propagation Equation

In our pump calibration example it's fairly obvious that the volumetric flow rate $$Y$$ (mL/s) can be calculated as the ratio of the measured volume $$V$$ (mL) to the discharge time $$t$$ (s), or $$Y = V/t$$. According to Eq. 1, the error of the volumetric flow rate is therefore

$e_Y^2 = \left(\frac{\partial Y}{\partial V}\right)^2 e_V^2 + \left(\frac{\partial Y}{\partial t}\right)^2 e_t^2,$

where $$e_V$$ and $$e_t$$ are the errors of the measured volume and measured time. Evaluating the partial derivatives yields

$e_Y^2 = \left(\frac{1}{t}\right)^2 e_V^2 + \left(-\frac{V}{t^2}\right)^2 e_t^2,$

or, after a bit of rearranging,

$\left(\frac{e_Y}{Y}\right)^2 = \left(\frac{e_V}{V}\right)^2 + \left(\frac{e_t}{t}\right)^2,$

where $$Y$$ is the calculated flow rate for the measured volume $$V$$ and measured time $$t$$.

### Example 2: Use of the Propagation Equation

For a knob setting of 1.0 +/- 0.2, the measured volume was 23 +\- 2 mL and the measured time was 10.0 +/- 0.2 s. Using Eq. 3 we calculate $$Y = 2.3 \pm 0.3$$ mL/s (don't forget to round up the individual errors before propagating them).

Application of the propagation equation derived in Example 1 to the data collected in Table 1, following the reporting rules summarized in Reporting the Error, and using the sources of error noted in (surprise) Sources of Error for the remaining values yields the data in Table 2.

Table 2. Propagation of error applied to experimental data for the hypothetical pump calibration experiment.
Knob Setting Discharge rate (mL/s)
0.0 +/- 0.2 0.0 +/- 0.0
1.0 +/- 0.2 2.3 +/- 0.3
2.0 +/- 0.2 4.2 +/- 0.4
3.0 +/- 0.2 7.2 +/- 0.5
4.0 +/- 0.2 8.4 +/- 0.6
5.0 +/- 0.2 11.9 +/- 0.7

The following MATLAB code was used to plot these data in Figure 1. You may wish to begin a MATLAB script with this material as I'll be using it in subsequent examples.

 knob = [0; 1; 2; 3; 4; 5];
rate = [0.0; 2.3; 4.2; 7.2; 8.4; 11.9];
err_knob = [0.2; 0.2; 0.2; 0.2; 0.2; 0.2];
err_rate = [0.0; 0.3; 0.4; 0.5; 0.6; 0.7];

figure; hold on;
h1 = errorbar(knob, rate, err_rate, err_rate, err_knob, err_knob, 'ko')
axis square
xlabel('Knob Setting');
ylabel('Discharge Rate (mL/s)');

## Regression Analysis

I assume that you've all taken your basic statistics course(s) and therefore I won't go into the details of linear and non-linear regression by least-squares analysis. Instead, I'll focus on how to use MATLAB to perform the most common statistical analyses that you'll need in this course.

To perform regression analyses in MATLAB there are two options:

1. The Curve Fitting App is an interactive feature of MATLAB that lets you play with various fitting procedures. It's quite helpful if you're not sure what form best represents your data or if you need to do some quick analyses and don't care about integration of analysis results into other functions.
2. The fit() function is the function underneath the hood of the Curve Fitting App and therefore has all the same functionality but you'll need to know how to use it and which functions it calls on to make various calculations. The examples below show you many of the common functions you'll need for this course and therefore I recommend you use fit() instead of the App for most analyses.

### Example 3: Linear Regression

The calibration plot in Figure 1 looks rather like a straight line so we'll start with a basic linear fit for our calibration curve. Assuming you've completed the code in the previous code block, the following code demonstrates how to perform a linear fit of our calibration data:

 >> fobj = fit(knob, rate, 'poly1')

fobj =

Linear model Poly1:
fobj(x) = p1*x + p2
Coefficients (with 95% confidence bounds):
p1 =  2.309   (1.976, 2.642)
p2 = -0.1048  (-1.113, 0.9035)

The fit object fobj is MATLAB object that holds all the information about our calibration curve. If, for example, we'd like to know the predicted flow rate for a knob setting of 2.5 then we simply pass 2.5 as an argument to the fit object,

>> fobj(2.5)

which gives a flow rate of 5.6667 mL/s. We need to compute the error on such a prediction--it's obviously not accurate to 4 digits of precision--but more on that below. To overlay our fitted equation atop the experimental data we create first a series of x-axis points for many arbitrary knob settings, calculate the corresponding predicted flow rate at each such value, and then overlay the plot. Adding to the previous MATLAB script we therefore have

 fknob = linspace(min(knob), max(knob), 50);
frate = fobj(fknob);
h2 = plot(fknob, frate, 'k-');

After these additions Figure 2 is produced.

### Example 4: Non-Linear Regression

Non-linear regression refers to a regression wherein the response variable depends in a non-linear way on the regression coefficients. A fit of the form $$y(x) = c_1 \textrm{exp}(c_2 x)$$ is non-linear because $$c_2$$ appears in the exponential but the form $$y(x) = c_1 \cos (x) + c_2 \sin (x)$$ is linear even though $$\sin (x)$$ and $$\cos (x)$$ are non-linear because $$y(x)$$ varies only linearly with the coefficients $$c_1$$ and $$c_2$$.

If for some reason we suspected the pump calibration curve varied in an exponential form then we can pass such a form to MATLAB's fit() function along with a starting guess for the coefficients (usually 1, 1 unless you have a better guess). The result is shown in Figure 3 below.

 >> fobj = fit(knob, rate, 'a*exp(b*x)', 'Start', [1 1])

fobj =

General model:
fobj(x) = a*exp(b*x)
Coefficients (with 95% confidence bounds):
a =  1.814   (4.033, 3.224)
b =  0.3833  (0.2041, 0.5625)

### Example 5: Custom Equations

Not all equations will match the functional form of the default equations, nor can they all be easily expressed by the simple character syntax used in Example 4 above. For more complicated equations we can instead define a custom function and tell MATLAB to use that instead.

Note 1: If you're comfortable with MATLAB then it's often faster to use anonymous functions for this process but I'm going to stick with the standard "separate function" system because it works in both cases and is somewhat more straightforward to understand.

Note 2: You can download the data set for this example here; it's not the same data set that the previous examples have used. Save it to your working directory and load it with load('pH_OLR_sampleData.mat').

As an example of a more complicated function consider the pH step-response $$Y(t)$$ of a first-order system with integrator[3],

$Y\left(t\right) = MK_p\left[t-\tau\left(1-e^{-t/\tau}\right)\right],$

where $$M$$ is the step size, and $$K_p$$ and $$\tau,$$ are coefficients to be determined by non-linear regression (these values will make more sense after you complete the Pre-Lab for pH Control, but for this level of detail is sufficient for our purposes here).

The first thing to do is create a function that returns a value for $$Y$$ given $$t, K_p$$ and $$\tau$$, and for this example we'll assume that $$M = 40$$. Note that the fit() function requires the dependent variable (here, $$Y(t)$$) to be defined as y and the independent variable to be assigned x.

 function Y = myFitFunc(x, K, tau)

Y = K .* 40 .* (x - tau .* (1 - exp(-x ./ tau)));

Save this function to your working directory, and now we use the fittype() function to tell MATLAB that myFitFunc() is something we want to use in a fitting procedure:

 >> ft = fittype('myFitFunc(x, K, tau)');

The data are assigned to more obvious variables to make the fit() call itself a little more intuitive and then we call fit() in a nearly analogous way as in Example 4 except with our custom fit type object instead of a custom equation.

 >> t = dat(:, 1); pH = dat(:, 2) - dat(1, 2);
>> Y = fit(t, pH, ft, 'StartPoint', [1 1])

Y =

General model:
Y(x) = myFitFunc(x, K, tau)
Coefficients (with 95% confidence bounds):
K =     0.02299  (0.01743, 0.02855)
tau =     4.784  (2.921, 6.647)

That's all there is to it! This is the most versatile way to call fit() because you can use it with almost any function you can think up, but of course most of the common ones are built-in and even the less common ones can frequently be expressed in the simple syntax of Example 4.

### Choosing the Best Form

Looking at Figure 3 the linear fit is obviously better than the exponential, but how can we quantify this recognition? We'll need stronger evidence than "it looks better," especially if the predictions are more similar than linear vs. exponential (for example, a higher-order polynomial will give results very similar to the linear fit).

#### Goodness-of-Fit

For a quick comparison the goodness-of-fit parameter, also known as $$R^2$$, is sufficient. More detailed and accurate evaluations can be obtained by ANOVA, specifically by comparing the P-value of the two models but for ease of use (i.e., to keep using fit()) we'll stick with the goodness-of-fit parameter. The $$R^2$$ parameter always varies between 0 and 1, and has the following interpretation:

• $$R^2 = 0$$ means that the model explains none of the variability of the response data around its mean.
• $$R^2 = 1$$ means that the model explains all of the variability around its mean.

To extract $$R^2$$ from the fit we add an output variable gof to our fitting procedure,

>> [fobj, gof] = fit(knob, rate, 'poly1');

and the $$R^2$$ value is returned as gof.rsquare. For the linear fit $$R^2 = 0.9893$$ and for the exponential fit $$R^2 = 0.9379$$, so we conclude that the linear fit is superior because more of the response data variability is explained by the linear model (specifically, 98.93% of the variability) than by the exponential model (93.79%).

#### Confidence Intervals

Using $$R^2$$ is fine if the results are different by more than a few percent as in the previous example, but what if two fits give comparable $$R^2$$ values? For example, a fit of the same data to a quadratic polynomial--achieved by using poly2 instead of poly1 in the call to fit()--yields $$R^2 = 0.9906$$, which is nearly the same as the linear model. When the models are so similar we instead rely on the confidence intervals of the coefficients to determine which are adding something meaningful and which can be discarded. The key principle is that if the coefficient's confidence interval contains 0 then the coefficient can be excluded.

The 95% confidence intervals are automatically reported as part of the fit calculation:

 >> fobj = fit(knob, rate, 'poly1')

fobj =

Linear model Poly1:
fobj(x) = p1*x + p2
Coefficients (with 95% confidence bounds):
p1 =  2.309   (1.976, 2.642)
p2 = -0.1048  (-1.113, 0.9035)

>> fobj = fit(knob, rate, 'poly2')

fobj =

Linear model Poly2:
fobj(x) = p1*x^2 + p2*x + p3
Coefficients (with 95% confidence bounds):
p1 =  0.05714  (-0.2258, 0.3401)
p2 =  2.023    (0.5492, 3.497)
p3 =  0.08571  (-1.481, 1.652)

There are two interesting features of these results:

1. The coefficient confidence intervals for p1 and p3 of the quadratic fit both contain 0 and therefore both p1 and p3 can be excluded from the fit. When you find that a coefficient can be excluded, don't simply ignore it and report the other coefficients; repeat the analysis for the new functional form so that the new form more accurately calculates the remaining coefficients. For example, since the quadratic form contains too many coefficients we'd reduce the order to a linear form.
2. The intercept (p2) confidence interval for the linear fit also contains zero, so we can simplify the model even further by eliminating the intercept.

For this particular calibration curve, excluding the intercept (i.e., setting it equal to zero) is acceptable because a 0 value on the pump corresponds to 0 flow rate. This may not be the case for all instruments, some of which may have a non-zero offset resulting in a non-zero intercept!

### Error Analysis

Generally you'll be using regression on a given data set to achieve one of two results:

1. You want the coefficients. This is the case when, for example, you're trying to fit the data to a known theoretical form and the coefficients are related to other physical quantities (e.g., reverse osmosis) or are themselves the values of interest (e.g., plate heat exchanger). In this case you need to determine the Error of Fitted Coefficients and you need to include the Error of Fitted Curves on your plots.
2. You want the prediction. This is the case when you're using the data as a calibration curve and you want to be able to use the resulting equation to predict new values of some variable (e.g., new flow rates for new pump knob settings). In this case you need to include the Error of Predicted Values on your plot.

Depending on what you're using the fit for, you'll need different error analyses as summarized below.

#### Error Of Fitted Coefficients

This is the easy one because these confidence intervals are provided as the 95% confidence bounds of a fit object. For example, if for some reason the slope of the pump calibration curve was the quantity of interest (and I doubt this would ever be the case) then you'd report the slope as $$b = 2.3 \pm 0.3$$, being careful to follow the reporting guidelines summarized above. You could then use this value and uncertainty in subsequent propagation-of-error formula if necessary.

#### Error Of Predicted Values

Any regression can also be used to produce prediction intervals, which are similar to confidence intervals but also account for the random variation expected in future data values, and are therefore bigger than confidence intervals. Such intervals are commonly used when the regression is used to produce a calibration curve, as we've done with the pump calibration example above, and you'd like to know the uncertainty of the predicted value (e.g., the uncertainty of the flow rate for a new pump knob setting).

Prediction intervals are calculated as

$y = \hat{y} \pm st_{\alpha, N-2} \left[1 + \frac{1}{N} + \frac{n\left(x - \bar{x}\right)^2}{n \sum x_i^2 - \left(\sum x_i\right)^2}\right]^{1/2},$

where $$\hat{y}$$ is the $$y$$ value predicted by the fit, $$t_{\alpha/2, N-2}$$ is the Student t value for a confidence level of $$\alpha$$ with $$N-2$$ degrees of freedom, $$N$$ is the number of data pairs, and $$s$$ is the residual standard deviation (also known as the residual standard error) available from ANOVA and calculated as

$s = \left[\frac{\sum\left(x_i - \bar{x}\right)^2}{N-2}\right]^{1/2}.$

Notice that in addition to being super ugly, the prediction interval--everything to the right of the $$\pm$$ symbol in Eq. (5)--depends on the $$x$$ value for which it is being calculated. Fortunately, MATLAB provides the predint() function to handle all of this nastiness with minimal fuss.

For example, let's say that you set the pump knob to 2.5. We saw earlier that the linear fit predicts a flow rate of 5.6667 mL/s, and now we'll use the prediction interval to calculate the uncertainty of this value:

>> predint(fobj, 2.5)

ans =

4.1619    7.1714

We would therefore report our best estimate of the flow rate as $$5.6667 \pm 1.5047 \textrm{ mL/s}$$, or more appropriately as $$6 \pm 2 \textrm{ mL/s}$$ since we need to follow the reporting guidelines.

Figure 4 illustrates several details of the prediction interval

• The prediction interval of a new measurement (in blue) is larger than the confidence intervals of the individual measurements which were used to create the calibration curve in the first place. This is because the new measurement has two sources of error: the random error of the value itself, as well as the error introduced by using a calibration curve.
• The prediction interval is curved with its narrowest region near the mean value of all calibration data points (usually near the center, if the data are uniformly spaced). This implies that your predictions will be less accurate if they fall near the extrema of your calibration curve, so plan accordingly!

The code to produce this figure is provided below.

 clear all
close all
knob = [0; 1; 2; 3; 4; 5];
rate = [0.0; 2.3; 4.2; 7.2; 8.4; 11.9];
err_knob = [0.2; 0.2; 0.2; 0.2; 0.2; 0.2];
err_rate = [0.0; 0.3; 0.4; 0.5; 0.6; 0.7];

figure; hold on;
h1 = errorbar(knob, rate, err_rate, err_rate, err_knob, err_knob, 'ko')
axis square
xlabel('Knob Setting');
ylabel('Discharge Rate (mL/s)');

fobj = fit(knob, rate, 'poly1');
fknob = linspace(min(knob), max(knob), 50);
frate = fobj(fknob);
h2 = plot(fknob, frate, 'k-');
fpi = predint(fobj, fknob);
h3 = plot(fknob, fpi, 'r-');
legend([h1 h2 h3(1)], {'Data', 'Fit', 'Pred. Int'}, 'Location', 'NorthWest');

k_new = 2.5;
r_new = fobj(k_new);
pi_meas = predint(fobj, k_new);
plot(k_new, r_new, 'bo', 'MarkerFaceColor', 'b');
plot([k_new, k_new], pi_meas, 'b-');

#### Error of Fitted Curves

As we saw above, if we fit a set of experimental data to a theoretical form then we expect a certain amount of uncertainty on the fitted coefficients; we called this uncertainty the "95% confidence intervals on the coefficients." For the purposes of visualization (e.g., on a plot) these numbers aren't particularly helpful and instead we use a different kind of confidence interval called a "functional confidence interval," sometimes also called the "confidence interval of the fit" or the "fit confidence interval."

The equation to determine the fit confidence interval at any point on the regression curve is

$y = \hat{y} \pm st_{\alpha, N-2} \left[\frac{1}{N} + \frac{n\left(x - \bar{x}\right)^2}{n \sum x_i^2 - \left(\sum x_i\right)^2}\right]^{1/2},$

which is identical to the equation for the prediction interval except the "1" inside the square root has been omitted. As with the prediction interval, MATLAB gives us a built-in function to calculate this value so that we don't have to deal with this messiness. Even better, it's the same predint() function we already know and love but with a few extra settings, namely predint(fobj, xvals, level, intopt, simopt), where fobj is the fit object, xvals are the x-values at which to calculate the confidence interval, level is the desired confidence level (always 0.95 for this course), intopt is an interval option that tells MATLAB which interval--prediction of confidence--to compute, and simopt is an option that tells MATLAB to calculate the individual or simultaneous confidence intervals (we always use simultaneous).

For example, to calculate the fit confidence interval of calibration data (something you should not do, but it's convenient to do so here since we've already got all the data) we would use the following code:

ci_fit = predint(fobj, fknob, 0.95, 'functional', 'on');
h4 = plot(fknob, ci_fit, 'b-');
legend([h1 h2 h3(1) h4(1)], {'Data', 'Fit', '95% Pred. Int.', '95% Conf. Int.'}, ...
'Location', 'NorthWest');
hold off

This code will produce the plot shown in Figure 5.

Note that the confidence interval is smaller than the prediction interval. This is because the confidence interval does not include the (predicted) random variation of a new variable; it is concerned only with the values already recorded. This is why I mentioned above that the fit confidence interval should not be included in calibration data, which are always concerned with the next data point. Fit confidence intervals, as described on the Plots and Tables page, should always be used when you perform a regression to compare experimental data to a theoretical form (e.g., rate laws in UV photocatalysis, characteristic curves in PEM fuel cells, mass transfer coefficients in reverse osmosis, and non-dimensional analyses in plate heat exchangers).

## Interpolation

Interpolation is the tool we use when we have a table of data but the value we seek isn't explicitly listed in the table. For example, consider this abbreviated enthalpy table for steam (in kJ/kg) shown in Table 3.

Table 3. Enthalpy table for steam.
P

(bar)

H (T = 100 °C)

(kJ/kg)

H (T = 150 °C)

(kJ/kg)

H (T = 200 °C)

(kJ/kg)

0.5 2683 2780 2878
1.0 2676 2776 2875
5.0 419 632 2855

Suppose we need to know the enthalpy of steam at 1.0 bar and 163 °C: there's a row for 1.0 bar but no column for 163 °C, so what do we do? We interpolate!

### Linear Interpolation

Linear interpolation is the most common form of interpolation and it assumes that the points given in a table are "connected" by straight lines. For example, consider the P = 1.0 bar row of Table 3 in graphical form in Figure 6.

The white circles are the data points given in Table 3, the red dashed line is a straight line between each point, and the red data point is the value at 163 °C that we wish to determine. Since the points are assumed to be connected by straight lines we can use similar triangles to determine the enthalpy value at 163 °C; graphically, the two triangles we need are shown in Figure 7.

The "similar triangles" principle applied to this situation is that the ratio of the arms of each triangle will be equivalent such that

$\left(\frac{\textrm{rise}}{\textrm{run}} \right )_{\textrm{triangle 1}} = \left(\frac{\textrm{rise}}{\textrm{run}} \right )_{\textrm{triangle 2}}$

Importantly, we know the x and y-values of all points in triangle 2 and all but one value for triangle 1, so we have one equation for one unknown:

$$\frac{2875-2776}{200-150}=\frac{H-2776}{163-150}$$

From this equation we calculate the unknown enthalpy as H = 2801.7 kJ/kg and this is what we call the interpolated value.

Another way to think of this process (and the one that I usually use) is to imagine a miniature table that contains only the points of interest:

T

(°C)

H

(kJ/kg)

150 2776
163 HT
200 2875

The "similar triangles" thing still applies here such that (163-150)/(200-150) has to be equal to (?-2776)/(2875-2776), from which we can again determine the enthalpy as HT = 2801.7 kJ/kg at 163 °C and 1 bar.

### Double Interpolation

If the data pair you seek isn't listed on either the rows or columns of the table then you'll need to perform double interpolation. For example, suppose you've got the enthalpy table for steam (in kJ/kg) given in Table 3 but you need the enthalpy for P = 0.77 bar and T = 179 °C. The data given in the table cover this range but this particular temperature and pressure are not given. You could draw the similar triangles thing again but I think it's easier to show in in the tabular form as

P

(bar)

H (T = 150 °C)

(kJ/kg)

H (T = 179 °C)

(kJ/kg)

H (T = 200 °C)

(kJ/kg)

0.5 2780 H3 2878
0.77 H1 HTP H2
1.0 2776 H4 2875

This table is kind of like Table 3 but "stretched out" a bit by the unknown enthalpies H1, H2, H3, and H4 between each of the tabulated values, as well as the central HTP which is the one we actually want. We have two options:

• Option 1: Interpolate to get H1, then again for H2, and finally a third time for HTP.
• Option 2: Interpolate to get H3, then again for H4, and finally a third time for HTP.

You should end up with the same answer regardless of which way you choose. For example, if we go with Option 1 then the table looks like

P

(bar)

H (T = 150 °C)

(kJ/kg)

H (T = 179 °C)

(kJ/kg)

H (T = 200 °C)

(kJ/kg)

0.5 2780 2878
0.77 2778.8 HTP 2876.4
1.0 2776 2875

The interpolation for HTP is (HTP - 2778.8) / (2876.4 - 2778.8) = (179 - 150) / (200 - 150) from which HTP = 2835 kJ/kg.

Option 2 yields the following:

P

(bar)

H (T = 150 °C)

(kJ/kg)

H (T = 179 °C)

(kJ/kg)

H (T = 200 °C)

(kJ/kg)

0.5 2780 2836.8 2878
0.77 HTP
1.0 2776 2833.4 2875

The interpolation for HTP is (HTP - 2836.8) / (2833.4 - 2836.8) = (0.77 - 0.5) / (1.0 - 0.5) from which HTP = 2835 kJ/kg. The reason for the equality is because the interpolation is linear: any point interpolated between the four tabulated points will always be on the plane defined by the four points, as illustrated in Figure 8.

### Using MATLAB for interpolation

Naturally MATLAB makes this whole process easier thanks to the interp1 and interp2 functions. You can use the documentation to check the syntax but if we want to repeat the single interpolation example we can do something like

T = [100 150 200];
H = [2676 2776 2875];
H_T = interp1(T, P, 163)

which yields the same value calculated above. To do the double interpolation is

T = [100 150 200];
P = [0.5 1.0 5.0];
H = [2683 2780 2878;
2676 2776 2875;
419  632 2855];
H_TP = interp2(T, P, H, 179, 0.77)

which again yields the same value calculated above. These are quite a bit faster than doing the interpolation by hand (I usually use Excel if I'm going to do that) but it can take a little while to get comfortable with the syntax and what the function is actually doing in order for this to be faster.

### Limitations

Sometimes it's not obvious whether you should use interpolation or regression to determine missing values, but a good rule of thumb is that if the data is of high precision and accuracy (like in the steam tables) then you should use interpolation instead of regression because you'll lose valuable information by using regression (try it with a steam table some time if you're bored).

Also, look again at the row for P = 5.0 bar in Table 3: why are the first two entries so much smaller than the other two? The answer is that these two data points correspond to a liquid phase whereas the other entries correspond to vapor phases. The variation is not linear between these points because of the abrupt change due to the enthalpy of vaporization, and it's up to you to identify such limitations when you use this technique.

## References

1. You can calculate the 95% confidence interval for repeated measurements in MATLAB using ci = tinv(0.975, length(z) - 2) .* std(z) ./ sqrt(length(z)), where z is a vector of measured values.
2. Finlayson, B.A.; Biegler, L.T., Mathematics. In Chemical Engineers' Handbook, Perry, R., Chilton, C., Eds. 8th ed.; McGraw Hill: New York, 2008, p 3-86.
3. Yes, this is the answer to one of your Pre-Lab questions for pH Control but you should derive it for yourself as an exercise.
Community content is available under CC-BY-SA unless otherwise noted.


text:
# All in the Corners

How to pick a random unit vector

As a quantum theorist, I spend a lot of time thinking about high-dimensional spaces. These are the playgrounds for quantum many-body systems, and they are vast. The technical name is a Hilbert space, and it’s the space of complex vectors with the additional structure of a way to put vectors together (called an inner product).

Hilbert space is big (see “The Curse of Dimensionality”), but the usable area for quantum theory is often much smaller. This means we are stuck finding corners of a high-dimensional space that describe physically-relevant phenomena.

A very important tool in my field is entropy. It can characterize the amount of entanglement between quantum states, which lets us talk about systems that produce behaviour far away from our usual ideas.

Recently, my supervisor and I were preparing homework for a class I help teach, and we wanted students to investigate the amount of entanglement1 present in a random quantum state of N qubits. This was a numerical exercise, and to complete it they needed to produce random unit vectors (because quantum states are normalized). This got me thinking about how to actually choose random unit vectors.

In this essay, we’ll explore the various ways you can choose random vectors, and why some methods won’t give you truly random unit vectors. To get there, we’ll need to think about high-dimensional spaces, why certain sampling procedures can lead to “clumping”, and why the requirement of being a unit vector changes the geometry of sampling we should be considering. To keep things from veering off into me telling you to visualize things in 35-dimensional space, we will build off our knowledge of two and three dimensions, and also note where things can go wrong.

By the end, we will see that the question of picking a random unit vector is more complicated than it seems at first glance.

## What is random?

When we use the word “random”, we need to be careful that we are talking about the same thing. That’s because random is always with respect to a probability distribution. Picking from a uniform distribution (each outcome is weighted equally) versus a normal distribution (outcomes near the mean are favoured, then taper off as you go further out) is a very different affair. You shouldn’t be surprised to get values near the endpoints of your interval in the uniform case, but you should be surprised if you’re picking from a normal distribution.

All this to say: random needs context.

We want to pick random unit vectors. As such, it’s probably worth defining exactly what we mean. A vector for our purposes will just be a list of numbers, with the number of entries corresponding to our desired dimension. For three dimensions, a vector could look like (a, b, c), where a, b, and c are real numbers. A unit vector just means that the norm of the vector is one. In particular, the norm we’re using is the regular L2 norm, which means we look at the sum of squares of the components, so our vector would have norm (a2+b2+c2)1/2. To transform this into a unit vector, we just need to divide each component by the vector norm.

Finally, we get to our lovely word: random. Here, I want random to mean that, out of the possible unit vectors, each one is equally likely to be chosen. In other words: I want to pick a unit vector uniformly at random.

With that out of the way, let’s dive into strategies to implement this. If you were to ask me, my first response would be a procedure like this:

1. For each coordinate of your vector, choose a random number in some interval [-a,a] (it doesn’t matter what a is).
2. Rescale the vector you get by its norm, producing a unit vector.
3. Repeat for the desired number of random unit vectors.

Simple, right?

The problem is that it doesn’t work.

Oh, you will get unit vectors. But they won’t be uniformly sampled. Instead, you will get “clumps” in certain regions of the possible unit vectors.

To check my claim, I’ll run some simulations that generate vectors at random, and then plot the results. To keep things visible, I’ll start by doing this in two dimensions. This will give us an easy in to the problem.

The following plot is what I get when I generate two lists: one with vectors whose components are picked uniformly from the interval [-1,1] (the red dots), and the other when I pick the vectors in the same way and then rescale them to be of unit length. (As a quick note, they the blue points aren’t found from the red ones. I generated two separate lists, but when you’re taking 5000 points, it shouldn’t be a big deal.)

Do you notice the clumping?

Okay, maybe it’s not so obvious. I’ve made the dots somewhat transparent so that overlap would show up as darker regions. If you look at the four “diagonal” directions (North-West, North-East, South-West, and South-East), you will notice the blue is more prominent. That’s because there’s clumping in those regions. Those areas are denser than the regions at the regular cardinal directions (North, East, South, West).

Why does this happen?

To get a handle on this question, think about what we’re sampling when we use the uniform probability distribution. Any given coordinate is being chosen in an interval [-1,1], so in total our vector is being sampled from a box whose side length is 1−(−1)=2. For the plot above, this box is just a square, and you can see that the red dots do indeed sample the box at a pretty uniform rate (there’s no clumping that I can see).

When we sample uniformly from the coordinates, the space we sample is a box.

But what are we really trying to do when we want to sample unit vectors?

Not only do we want to pick the coordinates at random, we want the norm of the vector to be one. This modifies the geometry of our space. If the vector must have unit length, the space we want to sample should be spherical. Looking at any other space will introduce bias in our results.

Okay, let’s unpack this a bit.

A unit vector is a vector with length one. This means that if your vector has coordinates xi, then the coordinates satisfy:

i xi2 = 1.

If we’re in two dimensions, this is just a fancy way of writing out the Pythagorean theorem with hypotenuse of length 1. If you vary the coordinates of your two-dimensional unit vector, you will trace out a circle. In higher dimensions, we trace out spheres, which is why we call such an equation spherically symmetric.

This means that choosing a random unit vector with uniform probability requires us to sample from a spherically symmetric distribution. Or at the very least, if you don’t use a spherical geometry, whatever you do to rescale things should counter the bias you introduced. We’ll just focus on making our distribution symmetric to start with, and see why picking uniform Cartesian coordinates doesn’t work.

Sampling our Cartesian coordinates as indicated above gives us a square in two dimensions, which we can plot relative to the circle (the spherically symmetric unit vectors).

Imagine you generate a bunch of random points within the square. This is what generating uniformly random Cartesian coordinates gives us (with the blue circle being the points that are unit vectors):

Now, we want to modify those dots so they all become unit vectors. How do we do this? We rescale! Each vector has a certain distance from the origin, and it’s likely that this distance isn’t one, like so.

Using the coordinates of the point, we can calculate it’s length from the origin, then divide each coordinate by that length. Doing so will give you a unit vector by definition.

But what does this rescaling look like on our square of points?

We know that rescaling will either push a vector with length less than one or pull a vector with length greater than one onto the circle.

Here is the fun part. Look at the difference between the North-East and North directions. In the North, almost all the points are closer to the origin than the circle. In fact, if you’re looking exactly North, than the only points that aren’t closer to the origin are those that already lie on the circle at the point (0,1). Everything else is pushed upwards.

But look at what happens in the North-East. You still have the same amount of points that are closer to the origin than the circle (because the circle is spherically symmetric, this is true in any direction you look). However, now you have a ton of extra points that are further from the origin than the circle!

I like to think of these as “extra” points that get pulled in because of the geometry of the box. In essence, the box puts “more” points in the corners, so those get rescaled down and populated the diagonal directions with more points when we rescale.

This is exactly what I was referring to when I talked about “clumping”. Because there are a bunch of extra points in the corners, your sampling will end up producing more points in those directions than the regular cardinal directions.

Here’s another way to think about. If you look in the North direction, you can go at most 1 unit away before you hit the box. This means that you have one unit of “space” to generate points on that line. But in the diagonal directions, you now have √2 units for your points to land on! Since your line is longer, more points will accumulate on it, making those directions favoured (biased).

When you use Cartesian coordinates that are sampled uniformly, you end up generating many more “corner samples” than you would want.

And it gets worse. How much “extra space” is there in a box compared to a sphere? In two dimensions, the area of a unit circle is π, and the box has area 22 = 4, so the ratio is π/4≈0.785. This means that there’s about 21.5% “extra space” between them.

In three dimensions, the volume of a unit sphere is (4/3)π, and the box’s volume is 23 = 8, so the ratio is (4/3)π/8≈0.524. As we can see, we have a lot more extra space!

If we look at the ratio in n dimensions, we need to know the sphere’s volume. Using Stirling’s formula (an approximation for factorials), the volume of a unit n-sphere is:

Vs(n) ~ (1/nπ)1/2(2πe/n)n/2.

We won’t dive into the reasoning here (perhaps in another essay). The volume of the box is much simpler. Each size as a factor of 2, so the volume is 2n. Therefore, the ratio looks like:

Vs(n)/2~ (1/nπ)1/2(πe/2n)n/2 → 0 as n → ∞.

Put more colorfully: As you go to higher dimensions, your box becomes almost entirely filled in the corners2. This means the bias we’ve seen in two dimensions will only get worse.

When you’re in high dimensions, everything is in the corners.

## The normal way

So what’s a good way to generate these unit vectors then?

One answer lies in taking a detour to see a certain kind of distribution, called a Gaussian or normal distribution. This distribution is one you see all over probability theory. In fact, I’d wager that saying “probability theory” will generate thoughts of the normal distribution in the minds of mathematicians.

What is the normal distribution? Briefly, the idea is that most of the possible values the distribution can take will be centered at a mean μ, and the the probability of getting a value away from the mean tapers off (controlled by a parameter σ called the variance). Mathematically, it looks like:

$$f(x) = \frac{1}{\sigma \sqrt{2\pi}} \exp{\left[ -\frac{1}{2} \left( \frac{x-\mu}{\sigma}\right)^2 \right]}.$$ This may look like an intimidating function, but when you plot it, the reality is much simpler.

The wonderful thing about a function like the exponential is that it transforms multiplication into addition. This is one of the laws that get taught to secondary students early on:

eaeb=ea+b, where a and b are numbers.

This is going to be hugely helpful to us.

We’ve seen graphically that a uniform distribution doesn’t give us a spherically symmetric distribution. What makes this distribution different? After all, we’re still going to sample each Cartesian coordinate independently, so what’s the big deal with the normal?

It all has to do with the joint probability distribution.

Let’s go back to our two-dimensional case. Here, we have our x and y coordinates that we want to sample. We can don’t really care about rescaling after, but it needs to be done in a way that everything stays spherically symmetric.

The joint probability distribution simply tells us how to link together multiple probability distributions. Because we’re looking at independent coordinates, our joint probability distribution will satisfy:

$$f(x,y) = f(x)f(y).$$ In the uniform case, we have f(x) = f(y) = 1/2, where the 2 comes in because it’s the length of the interval. This nets us a “square” distribution for the coordinates, which doesn’t respect spherical symmetry.

On the other hand, if we take a normal distribution with mean zero (μ=0) and variance one (σ=0), the joint distribution will look like this:

$$f(x,y) = f(x)f(y) = \frac{1}{2\pi} \exp{\left[-\frac{1}{2} \left( x^2 + y^2 \right) \right]}.$$ Look at the argument of the exponential. Does the expression x2+y2 look familiar?

This is precisely the radius r in polar coordinates! If we make the substitution, we end up with:

$$f(x,y) = f(r).$$ This is the reason that sampling from a normal distribution works. The structure of the exponential transforms a multiplication from independent distributions to an addition in the exponential. Furthermore, the fact that we have the square of the variable lets us employ the trick of going to polar coordinates.

Now that our distribution is a function of r, we’re guaranteed for it to be spherically symmetric. This means we’ve successfully constructed a way to sample unit vectors! All we have to do is the following.

1. For each Cartesian coordinate, sample a number from the probability distribution N(0,1) (the normal distribution with μ=0 and σ=1).
2. Rescale the vector you generate by its magnitude.

And that’s it! By simply changing the distribution we sample from, we can encode the fact that we want a spherically symmetric distribution.

——

Probability is a subtle mathematical subject, with the answer you get often hinging on the assumptions you baked into your approach (which are often implicit!). We saw a small example here, but this shows up everywhere in mathematics, and what it teaches me is that we need to be careful when tackling a problem. It’s worth spelling out your assumptions, checking to see if there are unintended consequences, and mending them if necessary.

This essay showed up that asking a simple question (How can I generate random unit vectors?) can lead us astray if we don’t question our assumptions. The probability distribution we used on the first try “fills up the corners” of our hypercube, bringing with it unwanted clumpy-ness. It also showed us that high-dimensional objects have properties that might run against our initial guess.

So the next time you start looking at hyperspheres and hypercubes, remember: the volume goes all to the corners.

## References

• For more on the n-sphere and its geometrical properties, see this Wikipedia article. There are plenty of interesting nuggets in higher dimensions, and many of them are collected here.
• This StackOverflow question is how I got started when I wanted to implement this in Python.
1. The measure we were interested in is called the entanglement entropy. It works like this. If you have a quantum state over a composite system (like many particles), you first divide it into two parts. Then, you can compute the reduced density matrix on one of those states, which is just a way to figure out the quantum state on one part by itself. Then, if you calculate the von Neumann entropy of that state, you get the entanglement entropy of the original quantum state. It also turns out that once you separate the composite system into two parts, it doesn’t matter which one you use to calculate the entropy. The result will be the same.

2. If you think it’s weird that the volume of the sphere takes up less and less space in the box as you increase the number of dimensions, welcome to the oddities of high-dimensional spaces!


text:
http://iet.metastore.ingenta.com
1887

## Peer-to-peer energy market between electric vehicles

• Author(s):
• DOI:

$16.00 (plus tax if applicable) ##### Buy Knowledge Pack 10 chapters for$120.00
(plus taxes if applicable)

IET members benefit from discounts to all IET publications and free access to E&T Magazine. If you are an IET member, log in to your account and the discounts will automatically be applied.

Recommend Title Publication to library

You must fill out fields marked with: *

Librarian details
Name:*
Email:*
Name:*
Email:*
Department:*
Why are you recommending this title?
Select reason:

ICT for Electric Vehicle Integration with the Smart Grid — Recommend this title to your library

## Thank you

Electric vehicles (EVs) are a major component of future electric grids, both for the increase in electricity demand and the flexibility they can add to the grid. Vehicle-to-grid and vehicle-to-building pilots have been tested and some have been approved by grid operators, but the EVs' possibilities shall be further enhanced. In previous works, the authors proposed a peer-to-peer energy market between EVs that largely reduced the expenses of their costly day-charging. This chapter further expands the model by taking into account the long-term effects of the market, which reduce the impact of the electric grid prices forecast on the market. The ratio between EVs that offer energy and those that demand energy is shown to be a good indicator for the market price forecast. Almost all energy demand occurs in pairs zone-time in which the number of offering EVs is more than five times the number of demanding EVs, for which the market price is very close to the electricity price at night.

Chapter Contents:

• List of abbreviations
• List of parameters and variables
• 7.1 Introduction
• 7.2 Activity-based model
• 7.3 Consumption model and drivers classification
• 7.3.1 Consumption model
• 7.3.2 Drivers classification
• 7.4 Intermediate charging process optimization
• 7.5.1 Determining the final price for the P2P trading system each TAZ and each time period
• 7.5.3 Proposed algorithm for solving the P2P market trading
• 7.6 Results of the P2P energy market
• 7.6.1 Individual analysis of vehicles from sets A and B
• 7.6.2 Electricity price analysis at TAZ level
• 7.7 Long-term peer-to-peer energy market
• 7.8 Conclusions
• References

Inspec keywords:

Preview this chapter:

Peer-to-peer energy market between electric vehicles, Page 1 of 2

| /docserver/preview/fulltext/books/tr/pbtr016e/PBTR016E_ch7-1.gif /docserver/preview/fulltext/books/tr/pbtr016e/PBTR016E_ch7-2.gif

### Related content

content/books/10.1049/pbtr016e_ch7
pub_keyword,iet_inspecKeyword,pub_concept
6
6
This is a required field