Monday, September 2, 2019

Kepler

Conic Section

 

Let's start with a circle
$$
r^2 =a^2
$$

If we look at this on a rectangular grid it decomposes like
$$
x^2 + y^2 =a^2
$$

or

$$
\left(\frac{x}{a}\right)^2+ \left(\frac {y}{a}\right)^2 = 1
$$

Ellipse

We can deform this into an ellipse by scaling the $x$ and $y$ axes independently

$$
\left(\frac{x}{a}\right)^2+ \left(\frac {y}{b}\right)^2 = 1
$$

With this form for the ellipse, the center and the major and minor axes are exposed clearly, but the focus is not. Let's look at another description of the ellipse, this time keeping the focus in focus. Begin again with the circle, $r=l$.  Now, we deform it into an ellipse like so:

\begin{equation} \label{eq1}
r=\frac{l}{1+\epsilon\cos\theta}
\end{equation}

If $\epsilon = 0$ this is a circle, and if $\epsilon$ is a small number less than $1$, this describes an ellipse centered at the focus.

The distance $r$ goes through extrema which are exposed easily if we invent a new parameter $a$, and rewrite as

\begin{equation} \label{eq2}
r=\frac{a(1+\epsilon)(1-\epsilon)}{1+\epsilon\cos\theta}
\end{equation}

When $\theta =0$, the denominator is $1+\epsilon$, and $r$ achieves a minimum

$$
r_<=a(1-\epsilon)
$$

When $\theta =\pi$, the denominator is $1-\epsilon$, and $r$ achieves a maximum

$$
r_>=a(1+\epsilon)
$$

These are the apsidal distances.  The semi-major axis, then, is

$$
a=\frac{r_> + r_<}{2}
$$

The focal distance is $a-r_<$, or

$$
\frac{r_> + r_<}{2} - r_< = \frac{r_> - r_<}{2} = \epsilon a \equiv c
$$

It's interesting to compare $a^2$ with $c^2$

$$
a^2 = \frac{1}{4} ({r_<}^2 + {r_>}^2 + 2 r_<r_>)
$$

$$
c^2 = \frac{1}{4} ({r_<}^2 + {r_>}^2 - 2r_<r_>)
$$

If we define a new parameter $b$ by
$$
b^2  \equiv  a^2 - c^2 = r_<r_>
$$

then we have $a^2 = b^2 + c^2$, a right triangle.

Apparently, $a$ is the arithmetic average of the apsidal distances, and $b$ is their geometric average.

The form \eqref{eq1} for the ellipse handles the limiting case of the circle, as we have seen, but also dials smoothly into the hyperbolic case, scooping up the parabolic case along the way.

\begin{align*}
\epsilon &= 0 \text{  circle}\\
\epsilon &>0 \text{  ellipse}\\
\epsilon &= 1 \text{  parabola}\\
\epsilon &>1 \text{  hyperbola}\\
\end{align*}

Hyperbola

One might wonder how the analysis that began with \eqref{eq2} would go through for the hyperbolic case. For that scenario we would, since $\epsilon >1$, invent a parameter $a$ and rewrite \eqref{eq1} like so

\begin{equation}
r=\frac{a(\epsilon+1)(\epsilon-1)}{1+\epsilon\cos\theta}
\end{equation}

Evidently, if $\theta = 0$, $r$ reaches a minimum of $r_<=a(\epsilon -1)$, while if $\cos\theta = -\frac{1}{\epsilon}$, $r$ goes bananas.  When $\theta$ reaches $\pi$, $r=r_>=a(\epsilon+1)$.  The midpoint between $r_<$ and $r_>$ is
$$
\frac{r_> - r_<}{2} = \frac {a (\epsilon +1) - a (\epsilon -1)}{2} = a
$$

while the focal distance is

$$
a + r_< =  \frac {a (\epsilon +1) + a (\epsilon -1)}{2} = a \epsilon \equiv c
$$

Similar to the ellipse case, if we define a new parameter $b$ by
$$
b^2  \equiv  c^2 -a^2 = r_<r_>
$$

then we have $a^2 + b^2 = c^2$, a right triangle.

Orbits

The motion of a mass in an inverse-square central force is described by a five-parameter family of curves.  You probably won't be surprised to learn that these curves are in fact the conic sections reviewed above.  After the plane (two parameters) and orientation (one parameter) are fixed, only two free parameters of the dynamics are available to describe the shape of the orbit, which will be some circle, ellipse, parabola, or hyperbola -- for instance the $a$ and $b$ used in the previous section.

A quick way to relate the parameters controlling the dynamics, to the parameters controlling the orbit is to consider the energy equation.

$$
E = \frac{m}{2} \dot{r}^2 + \frac{l^2}{2mr^2} - \frac{k}{r}
$$

Since the turning points of the motion describe the apsidal distances, we can set $\dot{r}$ to $0$ and read off the two apsides.

$$E =  \frac{l^2}{2mr^2} - \frac{k}{r}$$
$$Er^2 +kr -  \frac{l^2}{2m}  = 0$$

$$
r = -\frac{k}{2E} \pm \frac{\sqrt{k^2 + \frac{4El^2}{2m}}}{2E}
$$

$$
r = -\frac{k}{2E} \pm \frac{k}{2E}\sqrt{1 + \frac{2El^2}{mk^2}}
$$

Since the apsides are $a-c$ and $a+c$ (for the ellipse) we must have

$$ a= -\frac{k}{2E}$$
and
$$c= a\epsilon = -\frac{k}{2E}\sqrt{1 + \frac{2El^2}{mk^2}}$$

so that the major axis depends solely on the energy, and the focal distance depends on both the energy and the angular momentum.

Saturday, November 17, 2018

Entropy of a Conditional Distribution

Consider a set of $N$ objects, and let $g$ be the number of configurations indistinguishable under distribution $\rho$.  That is, the mulitplicity. Define $ \sigma$, the entropy, so that $g= e^{N\sigma}$. If $\rho$ is a multinomial distribution, entropy is asymptotic to

$\sigma = -\sum_i { p_i \ln{p_i}}$

If $\rho$ is a conditional distribution, $\rho(x|y)$, entropy looks like

$-\sum_i { p(x_i|y_j) \ln{p(x_i|y_j)}}$, which of course is a function of $y_j$.

The expectation value for this over all $y$ is
\begin{align}
\left< \sigma \right> &= -\left <\sum_i p(x_i|y_j) \ln{p(x_i|y_j)}\right > \\
 &= -\sum_j p(y_j) \sum_i  p(x_i|y_j) \ln{p(x_i|y_j)} \\
 &= \sum_j  \sum_i  p(y_j) p(x_i|y_j) [-\ln p(x_i, y_j) + \ln p(y_j)] \\
 &= \sum_j  \sum_i  (-p(x_i, y_j) \ln p(x_i, y_j) + p(x_i, y_j) \ln p(y_j)) \\
 &=  \sum_j  \sum_i -p(x_i, y_j) \ln p(x_i, y_j) + \sum_j p( y_j) \ln p(y_j) \\
&= \sigma(x,y) - \sigma(y)
\end{align}

In words: the conditional (on $y$) entropy is the joint entropy less the $y$-marginal entropy, or the conditional entropy plus the marginal entropy is the joint entropy.

$$
\sigma(x,y) =  \sigma(x|y) + \sigma(y)
$$

In a multivariate distribution, the sum of all the conditional entropies fails to sum to the joint entropy by a quantity known as the mutual information $I$

$$
\sigma(x,y) =  \sigma(y|x) + \sigma(x|y) + I(x,y)
$$

Entropy can be understood as the optimum average message length in the following sense. Consider a finite set of words, each with length $l_i$. Each word appears with probability $p_i$. The average word length, then, is $\sum_i p_i l_i$. Subject to constraint $\sum_i {2^{-l_i}} -1 =0$, the minimum length is achieved when $l_i = -\ln p_i$.

We can ask what the difference is between this optimal average length, and an average achieved with some other choice of $l_i$. Pick some non-optimal set of lengths $l_i = -\ln q_i$.  The difference, known as the Kullback-Leibler divergence, between the average of these lengths and the optimum average is

$$ \langle l \rangle - \sigma(\rho) = \sum_i p_i l_i - (- \sum_i p_i \ln p_i ) = \sum_i p_i \ln \frac{p_i}{q_i} $$

Mutual information can be viewed in terms of the KL divergence of the product distribution from the joint distribution like so

\begin{align}
D_{KL}(\rho(x,y)&||\rho(x)\rho(y)) = \sum_i \sum_j p(x_i, y_j) \ln \frac{p(x_i, y_j)}{p(x_i) \rho(y_j)} \\
&=  \sum_i \sum_j p(x_i, y_j) \ln p(x_i, y_j) -  \sum_i  p(x_i)  \ln p(x_i) - \sum_j p(y_j)  \ln p(y_j) \\
&= -\sigma(x,y) + \sigma(x) + \sigma(y) \\
&= -\frac12 (\sigma(x|y) + \sigma(y|x) + \sigma(x) + \sigma(y) ) + \sigma(x) +\sigma(y) \\
&= -\frac12 (\sigma(x|y) + \sigma(y|x) - \sigma(x) - \sigma(y) )  \\
&=  -\frac12 (-2I) = I
\end{align}

Monday, October 9, 2017

Representation Theory

Groups on Linear Spaces

In physical mathematics we frequently speak of group operations on linear spaces.  Consider the elements of a group ${G_a}$ and linear operators $T(G_a)$, so that $T(G_a)v \implies v'$, where $v$ and $v'$ are vectors.

\begin{align}
T(G_a) e_i &= T_{ji}(G_a) e_j \\
(e_j, T(G_a) e_i) &= T_{ji}(G_a) \\
T(G_a) T(G_b) &= T(G_c) \text{, if $ G_a G_b = G_c$}
\end{align}

The operators $T$ are an example of a representation of the group ${G_a}$ with respect to operator multiplication.

What is a group?

We restrict attention to finite groups, which are sets with a binary operation that is 
$\require{enclose}$
1.  closed, and
2.  associative

Also,
3.  $\exists$ an inverse for every element, and
4.  $\exists$ an identity element

Reducibility

This is an important topic in rep theory. Suppose L is a vector space, invariant with respect to the transformations $T(G_a)$. If $L=L_1 + L_2$, two orthogonal spaces also invariant with respect to $T(G_a)$, then T reduces. Otherwise, $T$ is irreducible. Amazing fact: all finite groups have a finite number of (non-equivalent) irreducible representations. We will show this.

Schur's Lemma

Before going further, we need a preliminary result.

$$
T^{\alpha}(G_a) A - A T^{\beta}(G_a) = 0 \; \enclose{horizontalstrike}{\lor}  \;G_a \implies A = \lambda \delta_{\alpha \beta} I
$$

But we can construct just such an $A$

$$
A = \sum_b T^{\alpha} (G_a) X T^{\beta}(G_a^{-1})
$$
$X$ is arbitrary, but determines $\lambda$

$$
\sum_{bmn} T_{im}^{\alpha}(G_b) X_{mn} T_{nj}^{\beta}(G_b^{-1}) = \lambda \delta_{\alpha \beta} \delta_{ij}
$$

Choose $X_{mn} = \delta_{mq}\delta_{np}$

\begin{equation}
\sum_{b} T_{iq}^{\alpha}(G_b)T_{pj}^{\beta}(G_b^{-1}) = \lambda \delta_{\alpha \beta} \delta_{ij} \label{eq:1}
\end{equation}

Let $\alpha = \beta$ and $i=j$

$$
\sum_{b} T_{iq}^{\alpha}(G_b)T_{pi}^{\beta}(G_b^{-1}) = \lambda
$$

Sum both sides over $i$

$$
\sum_{b} T_{pq}^{\alpha}(E) = \lambda S_{\alpha}
$$

$ S_{\alpha}$ is the dimension of the $\alpha$ representation.

\begin{equation}
g \delta_{pq} = \lambda S_{\alpha}  \label{eq:2}
\end{equation}


$g$ is the number of elements in the group.

Using $\eqref{eq:2}$ in $\eqref{eq:1}$:

\begin{equation}
\sum_{b} T_{iq}^{\alpha}(G_b)T_{pj}^{\beta}(G_b^{-1}) = \delta_{pq} \delta_{\alpha \beta} \delta_{ij} \frac{g}{S_{\alpha}} \label{eq:3}
\end{equation}
Let $T$ be unitary, $i=j$, $\alpha = \beta$,  and $p = q$.

\[
\sum_{b} T_{iq}^{\alpha}(G_b)T_{iq}^{* \alpha}(G_b) = \sum_{b}{|T_{iq}^{\alpha}(G_b)|}^2 = \frac{g}{S_{\alpha}}
\]

Characters

In $\eqref{eq:3}$, let $q=i$ and $p=j$ and sum over $i$ and $j$

\begin{align*}
\sum_{ijb} T_{ii}^{\alpha}(G_b)  T_{jj}^{* \beta}(G_b) & = \sum_{ij} \delta_{ij} \frac{g}{S_{\alpha}} \delta_{\alpha \beta} \\
\sum_{b} \chi^{\alpha}(G_b)  \chi^{* \beta}(G_b)  &= g \delta_{\alpha \beta}
\end{align*}

Where the character $\chi^{\alpha}(G_a)$ is the trace of the operator $T^{\alpha}(G_a)$

If we let $C_p$ count the members of a given class, the elements of which share a character, we can write this as

\begin{align}
\sum_p C_p \chi_p^{\alpha} \chi_p^{* \beta} &= g \delta_{\alpha \beta}  \label{eq:4} \\
\sum_p C_p |\chi_p|^2  &= g
\end{align}

For any representation, $\chi_p = \sum_i m_i \chi_p^i$ since any rep can be resolved to its irreducible components.  We can find the $m_i$, given the characters.

\begin{align*}
\sum_p C_p \chi_p^{\alpha}  \chi_p &= \sum_i \sum_p C_p \chi_p^{\alpha} m_i \chi_p^i \\
&= \sum_i g \delta_{\alpha i} m_i \\
&= m_\alpha g
\end{align*}

For any rep, $\sum_p |\chi_p|^2 \gt g$ with equality only if the rep is irreducible.

\begin{align}
\sum_p |\chi_p|^2 &= \sum_p C_p \sum_{\alpha} m_\alpha \chi_p^\alpha \sum_{\beta} m_\beta \chi_p^\beta \\
&= \sum_{\alpha \beta} \sum_p C_p  \chi_p^\alpha \chi_p^{* \beta}  m_\alpha  m_\beta \\
&= \sum_{\alpha \beta} \sum_p  g \delta_{\alpha \beta} m_\alpha  m_\beta \\
&= \sum_\alpha g |m_\alpha|^2
\end{align}

This equals $g$ only if the $m_\alpha$ are all $0$, with one exception, which is $1$

The Regular Representation


Define the $g$-dimensional representation $T^\text{Regular}$ by

$$
G_a G_b = \sum_c T_{cb}^\text{Regular}(G_a) G_c
$$

The $\chi^R(G)$ are all zero except $\chi^R(E)$ which is $g$.  Consider the reduction of the regular representation:

\begin{align}
\sum_p C_p \chi_p^{* \alpha} \chi_p^R &= g m_\alpha \\
C_E \chi_E^{* \alpha} \chi_E^R(E) &= g m_\alpha \\
S_\alpha g  &= g m_\alpha
\end{align}

So $m_\alpha = S_\alpha$, and $g=\sum_\alpha m_\alpha S_\alpha = \sum_\alpha S_\alpha^2$

The number of orthogonal vectors, $T_{ij}^\alpha$, is $\sum_\alpha S_\alpha^2$, which we now know is $g$. That is, these vectors span the vector space. We could expand any vector in this space as

$$
v = \sum_{ij\alpha}  Z(ij\alpha) T_{ij}^\alpha
$$

Or in component form
$$
v_a = \sum_{ij\alpha}  Z(ij\alpha) T_{ij}^\alpha(G_a)
$$

Character Relations

Now, consider a vector $v$ that has the same component along all directions in a class, like so:
$$
v_a = \frac {1}{g} \sum_{b=1}^g v_c  \,\,\, \text{ where $G_c = G_b^{-1} G_a G_b$}
$$

and decompose $v_a$ and $v_c$ along $T$ like so:

\begin{align}
v_c &= \sum_{ij\alpha}  Z(ij\alpha) T_{ij}^\alpha(G_c) \\v_a &= \sum_{b}  \sum_{ij \alpha} \sum_{kl} \frac{Z(ij\alpha)}{g} T_{ik}^\alpha(G_b^{-1}) T_{kl}^\alpha(G_a) T_{lj}^\alpha(G_b)
\end{align}

Using $\eqref{eq:3}$ to sum out $b$

\begin{align}
v_a &=    \sum_{ij \alpha}  \sum_{kl} \frac {Z(ij\alpha)}{g} \frac{g}{S_\alpha} \delta_{kl}  \delta_{ij}  T_{kl}^{\alpha}(G_a) \\
&= \sum_{i\alpha} \sum_k \frac{Z(ii\alpha)}{S_\alpha} T_{kk}^{\alpha}(G_a) \\
&= \sum_{i\alpha} \frac{Z(ii\alpha)}{S_\alpha} \chi^{\alpha}(G_a) = \sum_{\alpha} \chi^{\alpha}(G_a)  \left[ \sum_i \frac{Z(ii\alpha)} {S_\alpha} \right] = \sum_\alpha z(\alpha) \chi^\alpha(G_a)
\end{align}

We see that the $\chi^\alpha$ span this space of classes so we conclude there must be $n$ $\alpha$'s if there are $n$ classes.

\[
\bbox[5px,border:2px solid red]
{
\text{# of irreps = # of classes}
}
\]

Since the character table is square, we can re-write $\eqref{eq:4}$, which tells us the rows of the table are orthogonal, in the following form

$$
B_{\alpha p} = \sqrt{\frac{C_p}{g}} \chi_p^\alpha
$$

$$
\sum_p  B_{\alpha p} B_{\beta p}^{*} = \delta_{\alpha \beta}
$$


So $Det(B) = Det(B^{\dagger}) = 1$, so $B$ has an inverse, $B^{-1} = B^{\dagger}$, which implies the columns are also orthogonal:

$$
\sum_{\alpha} B_{\alpha p}^{*} B_{\alpha q} = \delta_{pq}
$$

$$
\sum_{\alpha} \frac{C_p }{g} \chi_p^{\alpha *} \chi_q^{\alpha} = \delta_{pq}

$$

\begin{equation}
\sum_\alpha \chi_p^\alpha \chi_q^\alpha = \frac{g}{C_p} \delta_{pq}
\end{equation}

Woo!

Projection Operators

$$
G_a e_i^{\alpha} = \sum_j T_{ji}^\alpha(G_a) e_j^{\alpha} \, \text {  (defines a representation) }
$$

\begin{align}
\sum_a  T_{mn}^{\beta}(G_a^{-1}) G_a e_i^\alpha  &= \sum_{aj} T_{mn}^{\beta}(G_a^{-1}) T_{ji}^\alpha(G_a) e_j^\alpha \\

&= \sum_j \frac{g}{S_\beta} \delta_{mi} \delta_{nj} \delta_{\alpha \beta} e_j^{\alpha} \\
&= \frac{g}{S_\beta} \delta_{mi}  \delta_{\alpha \beta} e_n^{\beta}
\end{align}

Let $m=n$
$$
\sum_a T_{nn}^{\beta}(G_a^{-1}) G_a e_i^\alpha = \frac{g}{S_\beta} \delta_{ni}  \delta_{\alpha \beta} e_n^{\beta}
$$

Define a projection operator

$$
\left[ \sum_a \frac{S_\beta}{g} T_{nn}^{\beta}(G_a^{-1}) G_a  \right] e_i^\alpha= \delta_{ni}  \delta_{\alpha \beta} e_n^{\beta}
$$

$$
P_n^\beta e_i^\alpha = \delta_{ni} \delta_{\alpha \beta} e_n^{\beta}

$$

Thursday, August 25, 2016

Bayes Thinking

Preamble

Consider a set $\{x_i\}$, and map $f(x)$ of that set: $y_i = f(x_i)$. Further, consider a set of measurements $\{d_i\}$ of that mapping. The measurements possibly differ from the mapped values $\{y_i\}$ like so:  $ d_i - \epsilon_i = y_i $. We don't know much about $\{\epsilon_i \}$, just that they are randomly selected from a gaussian distribution with center $0$ and width $\sigma^2$. So for a given $f(x_i)$, we see $d_i$ with probability $Normal(\epsilon_i)$.

Perhaps this is labeled data from a supervised learning scenario, or a regression problem. We wish to know which function, from the space of all possibilities, gives the best mapping. Bayes Theorem helps us understand how to think about a problem like this.

Bayes Theorem

\[ P(h \vert d) = \left( \frac{P(d \vert h)}{P(d)} \right) P(h) \]

Bayes theorem relates the probability of hypothesis $h$ given data $d$ to the probability of data $d$ given hypothesis $h$. The hypotheses are our candidate functions introduced above. These terms are sometimes described as $posterior = \left( {likelihood \over evidence} \right) {prior}$.  This is the math that tells you how to update your understanding from prior to posterior, given information likelihood, and  background evidence. Take a moment to appreciate what a bold statement that is. This tells you, essentially, how to think. Or at least how to learn.

So back to our problem: which is the best function? Clearly we want to choose the $h$ such that the probability for the winning $h$ is the highest in the sense of $ P(h \vert d)$. This is called maximum a-posteriori (MAP). We might write this as

\[ argmax \left\{ \left( \frac{P(d \vert h)}{P(d)} \right) P(h) \right\} \]

The normalizing factor $P(d)$ does not distinguish choice of hypothesis, and the same can be said of the prior if the prior is uniform.  With both of these assumptions in place, maximum a-posteriori is the same as maximum likelihood (ML):

\[ argmax \; P(h \vert d)  = argmax \; P(d \vert h) \]

Bayes Learning

So let's find

\begin{align*}
 h_{ML} &= argmax \; P(d \vert h)  \\
 &=   argmax \prod_i P(d_i \vert h) \\
 &=   argmax \prod_i Normal(\epsilon_i) \\
 &=   argmax \prod_i C \, \exp \left(-\frac{\epsilon_i^2}{2 \sigma^2} \right)
\end{align*}

We can take the log of the argmax expression, as the extremum is maintained.
\[ \bbox[5px,border:2px solid red] {
 h_{ML} =   argmin \sum_i \epsilon_i^2
} \tag{minimum squared error} \]

Whoever said that we should minimize squared error knew what he was talking about. Bayes agrees.

Let's back up a little and think about this another way. Consider the MAP hypothesis

\[ h_{MAP} = argmax \left\{ \left( \frac{P(d \vert h)}{P(d)} \right) P(h) \right\} \]

Again, $P(d)$ does not distinguish hypotheses, so we can drop that term

\[ h_{MAP} = argmax \;  P(d \vert h)  P(h)   \]

The maximum is maintained through the log function

\[ h_{MAP} = argmax [  \log P(d \vert h) + \log P(h) ]  \]

\[ h_{MAP} = argmin [  -\log P(d \vert h) -\log P(h) ]  \tag{Occam's razor}\]

With this relation in view, consider the first term. We decided earlier that was related to our measurement or calculation error (or both). The other term, $-\log P(h)$ can be considered a description (code) length, if we're thinking like C. Shannon. So we can paraphrase:

\[
\bbox[5px,border:2px solid red]
{
 h_{MAP} = \text{minimum error} + \text{minimum description length}
}
\]

Or, the best hypothsis is the simplest one that gets the job done.  Sounds like Occam's advice.

Wednesday, November 25, 2015

Some Results in Geometric Algebra

$\newcommand{\abs}[1] {\left\vert #1 \right\vert}$

Cramer's Rule

The roots of Cramer's rule have always seemed a little mysterious. If we have a system of equations

\[ a_1 x_1  + b_1 x_2 = c_1 \]
\[ a_2 x_1 + b_2 x_2 = c_2 \]

Then Cramer says
\[ x_1 = \frac{\det  \begin{bmatrix} c_1 & b_1 \\  c_2 & b_2  \end{bmatrix}}{\det  \begin{bmatrix} a_1 & b_1 \\  a_2 & b_2  \end{bmatrix}} \]

and

\[ x_2 = \frac{\det  \begin{bmatrix} a_1 & c_1 \\  a_2 & c_2  \end{bmatrix}}{\det  \begin{bmatrix} a_1 & b_1 \\  a_2 & b_2  \end{bmatrix}} \]

That determinants play any role in isolating linear parameters is not obvious. So it's delightful how geometric algebra shines light on this.

A note on notation: in this section, latin letters are vectors and greek letters are scalars, so that

\[ x = \alpha a + \beta b\]

is a vector equation relating vector $x$ to vectors $a$ and $b$. The gift we are given by geometric algebra is that vectors become invertible, so we can write, since $ a \wedge a = 0 $:

\[ a \wedge x = \beta a \wedge b \]
\[ \beta = \frac {a \wedge x}{a \wedge b} \]

And similarly,

\[ \alpha = \frac {x \wedge b}{a \wedge b} \]

so that

\[ x = \left( \frac {x \wedge b}{a \wedge b} \right) a +  \left( \frac {a \wedge x}{a \wedge b} \right) b \]

Now, in $\mathbb{R_2}$, 2-vectors are pseudoscalars, and pseudoscalars are scalar multiples of $I$, the unique (up to a scalar) highest-grade multivector. The determinant is just that scalar (say, $\lambda$): $a \wedge b = \lambda I$. So we can write:

\[ x = \left(\frac {\det (x \wedge b) I}{\det (a \wedge b)  I} \right) a +   \left( \frac {\det (a \wedge x) I}{\det (a \wedge b) I} \right) b \]


\[ x = \frac {\det (x \wedge b) }{\det (a \wedge b)  } a +   \frac {\det (a \wedge x) }{\det (a \wedge b) }  b \]

That's the fastest path to Cramer's Rule anywhere.

BAC-CAB

Consider the product of a vector and bivector: $ a (b \wedge c)$.

\[ \frac{1}{2} a (bc - cb) \]

\[ \frac{1}{2} (abc - acb) \]

Since $ab =2 a \cdot b -ba$, we can write

\[ \frac{1}{2} (2 a \cdot b \; c - bac -2a \cdot c \; b -cab ) \]

\begin{equation} \label{eq1}
a \cdot b c - a \cdot c b -\frac{1}{2} (bac-cab)
\end{equation}

The last two terms can be expanded as before:

\[ - \frac12(bac-cab) = -b a \cdot c + \frac{bca}{2}  + c a \cdot b - \frac{cba}{2}  \]

So in the expanded \eqref{eq1}, the dot products combine to give us

\[ 2  a \cdot b  c - 2  a \cdot c  b  + \frac{(bc-cb) a}{2} \]

So finally

\[ a (b \wedge c)  = 2  ( a \cdot b  c -  a \cdot c b ) + (b \wedge c) a \]

\[ \frac{ a (b \wedge c) - (b \wedge c) a}{2} =  a \cdot b c - a \cdot c b \]

\[ a \cdot (b \wedge c) = a \cdot b c - a \cdot c b \]

This is a handy result.  And, not restricted to $\mathbb{R^3}$.

Using the Dual to Exchange Dot and Wedge Products

We will develop a nice result with the dual operation that lets us exchange dot and wedge products.

For any multivector $A$, we can write $A = A^\perp I$, where $A^\bot$ is whatever multivector you need to multiply with $I$ to get $A$, called the dual of $A$. In other words, taking the dual is multiplying on right with $I^{-1}$. Now consider the dual of $AB$

\[ (AB)^\bot = AB^\bot\]

This can be confirmed by examining $(AB)I^{-1} = ABI^{-1}$.

Now compare grades of this expression, and conclude that the largest grades give us

\[ (A \cdot B)^{\bot} = A \wedge B^\bot \]

and the smallest grades give us

\[ (A \wedge B)^{\bot} = A \cdot B^\bot \]

Gibbs' Vector Product

The dual lets us move between Gibbs-style cross products $a \times b$ and geometric products:
$$  a \times b \equiv (a \wedge b)^{\bot} $$

Of course Gibbs' cross product is only defined in $ \mathbb{R^3}$.

Here's an example, which for some readers may complete the BAC-CAB story:
\begin{align}
a \times (b \times c) & = \{ a \wedge  (b \times c) \} ^\bot \\
& = a \cdot (b \times c)^\bot \\
& = -a \cdot b \wedge c
\end{align}
 
Another:
\begin{align}
a \cdot b \times c & = a \cdot (b\wedge c)^\bot \\
& = (a \wedge b \wedge c)^\bot
\end{align}

Rotations

On our way to describing angular momentum under a geometric algebra framework, we need to establish a preliminary object known as a rotor. And on the way to the rotor, we need to work through reflections. Consider a reflection axis $n$ ($n$ is a unit vector) and a vector $v$.

$$ v = v n n^{-1} $$
$$ = [ v \cdot n + v \wedge n ] n^{-1} $$

The first term is the component of $v$ along $n$.  So, the reflected $v$ is just

$$ v' =  -v \cdot n n^{-1} + v \wedge n  n^{-1} $$

A small manipulation lets us write

$$ v' =  -n \cdot v n^{-1} - n \wedge v  n^{-1} $$
$$  =  - [nv] n^{-1} = - nvn^{-1}  $$

Kinda cute. And, due to Hamilton, this gets us halfway to rotations, because he noticed that a rotation is simply two reflections. Say we reflect along $b$, followed by $a$:

$$ v_{rotated} = (ab)v(ab)^{-1} = RvR^{-1} $$

$R$ is known as a rotor: a product of two unit vectors.

$$ R = ab = a \cdot b + a \wedge b $$

We can find the magnitude of $ a \wedge b $ by writing

$$ ba =  b \cdot a + b \wedge a = a \cdot b - a \wedge b$$

and then considering the following product:

$$ abba = a^2 b^2 = 1 = (a \cdot b)^2 - (a \wedge b)^2 $$

Since $ a \cdot b $ is $\cos \theta$ we can write

$$ 1 - \cos^2 \theta = \sin^2 \theta = - (a \wedge b)^2 $$

So $ a \wedge b $ is a bivector whose magnitude is $\sin \theta $ and square is negative. Let $I$ be the unit bivector

$$ I = \frac{a \wedge b}{\sin \theta} $$

and write

$$ ab = R = \cos \theta + I \sin \theta = \exp (I \theta) $$

Looks like the common exponential phase factor $ e^{i \theta}$. Aside: what are the complex numbers?  In two dimensions, vector $(a,b)$ is related to complex $a + ib$ like so: $a \hat{x} + b \hat{y} = \hat{x} (a + \hat{x} \hat{y} b) = \hat{x} (a + I b)$. Just pre-multiply by $\hat{x}$, selecting the real axis. You always knew there was something like this going on behind the scenes.

Moving along: let's use the dagger operator to reverse the order of vector products: $ (ab)^\dagger = ba $. Since unit vectors are their own inverse, we can write $y$ as a rotated version of $x$ like so:

$$ y = R x R^{-1} = R x R^\dagger  $$

If $R$ is in fact a function of time, we can consider $ \dot{y} $:

\begin{align}
\dot{y} &= \dot{R} x R^\dagger + R x \dot{R^\dagger} \\
&= \dot{R} R^\dagger y + y R  \dot{R^\dagger} \\
&= (\dot{R} R^\dagger) y - y (\dot{R}R^\dagger)
\end{align}

For the last step, take the time derivative of $R R^\dagger = 1$ and show $\dot{R} R^\dagger = - R \dot{R}^\dagger$.

Let $ \Omega= - {1 \over 2} \dot{R} R^\dagger  $. $\Omega$ is apparently of even grade, and has no scalar part. $\Omega$ must be a bivector, and we can write:

$$ \dot{y} = - \Omega\cdot y = y \cdot \Omega$$

Orbital Dynamics

With rotational motion in our toolkit, let's develop the first few results in central-force dynamics. A nice trick for describing this is to shift the dynamics to a rotor. We accomplish this by considering a vector $x$ as a rotation of some basis vector, plus a stretch.

$$ x= U e_1 U^{\dagger} $$

We are now allowing a stretch as well as rotation, so U is a little more than a rotor: we are now considering a general even element. In three dimensions this is a scalar plus bivector [side note: the spinor and quaternion concepts fall out of this]. This presents us with a slight issue right away, in that U has four degrees of freedom, but we only need three for $x$. This allows us to apply a constraint on $U.$ But before deciding on that constraint, consider the rate of change of $x$:

$$ \dot{x} = \dot{U} e_1 U^\dagger + U e_1 \dot{U}^\dagger $$

Under the special case where $\dot{U} e_1 U^\dagger = U e_1 \dot{U}^\dagger $ (this is our constraint) we can write

$$ \dot{x} = 2 \dot{U} e_1 U^\dagger $$

Introduce a new variable s, such that $ {dt \over ds} =r= U U^\dagger$

$$ 2 {dU \over ds} = \dot{x} U e_1 $$
$$ 2 {d^2 U \over {ds}^2} = r \ddot{x} U e_1+ \dot{x} {dU \over ds} e_1 =  (\ddot{x} x+{1 \over 2} \dot{x}^2) U $$

The last step follows because $r  U e_1 = x U$. If we have inverse square central force motion

$$ m \ddot{x} = -{k x\over r^3} $$

then we have

$$ 2 {d^2 U \over {ds}^2} =  (\ddot{x} x+{1 \over 2} \dot{x}^2) U $$
$$= ( -{k \over m} {x^2 \over r^3} + {1 \over 2} \dot{x}^2) U $$
$$=  {1 \over m} ( -{k \over r} + {1 \over 2}m \dot{x}^2) U $$
$$ {d^2 U \over {ds}^2} = \left( {E \over 2m}  \right) U $$

Harmonic motion. Nice.

Orbit Shape

A mass $m$, at $r$, moving with respect to some origin, defines a plane described by the bivector $L$, the angular momentum:

\begin{align*}
L &\equiv r \wedge p \\
&= m r \wedge \dot{r}\\
&= m r \wedge (r \dot{\hat{r}} + \dot{r}\hat{r})\\
&= m r^2 (\hat{r} \wedge \dot{\hat{r}}) = m r^2 \hat{r} \dot{\hat{r}}
\end{align*}

For motion in a potential $V= -{k \over r}$ we have

$$ \dot{v} = -{k \over m r^2} \hat{r} $$

Form the product of this with $L$:

$$ L\dot{v} = -mr^2\dot{\hat{r}}\hat{r}  \left( {-k \over mr^2} \right) \hat{r} = k\dot{\hat{r}} $$

$$ {d \over {dt}} (Lv - k \hat{r}) = 0 $$

And we see that $Lv- k \hat{r}$ is a constant of the motion. This vector, which sits in the plane of motion $(L \wedge v = 0)$, is known as the Runge-Lenz vector. Let's be explicit about the conserved vector and call it $e$:
$$ Lv = k \hat{r} + ke$$

With this definition, $e$ is dimensionless. If we multiply through with $r$ we have

\begin{equation} \label{eq2}
Lvr = k \hat{r} r + ker 
\end{equation}

The left hand side is

\begin{align}
Lvr &= L (v \cdot r + v \wedge r)\\
&= L \left(v \cdot r + \frac{L^\dagger}{m} \right) \\
&= L (v \cdot r) + \frac{l^2}{m} \label{eq3}
\end{align}

The right hand side is

\begin{align*}
k \hat{r} r + ker &= k \abs{r} + k(e \cdot r + e \wedge r)   \\
&= k \abs{r} + \abs{e} \abs{r} \cos \theta + e \wedge r
\end{align*}

So we have, finally

$$ l^2/m + (r \cdot v) L = k \abs{r} + k \abs{e} \abs{r} \cos \theta + e \wedge r $$

The scalar part is

$$l^2/m = k \abs{r} + k \abs{e} \abs{r} \cos \theta $$

or

$$ \abs{r} = \frac{l^2/mk}{1+\abs{e} \cos \theta}$$

and the bivector part is

$$ (r \cdot v) L = e \wedge r  $$

or

$$ L = \frac{e \wedge r }{r \cdot v} $$

Returning to $\eqref{eq2}$ we can develop a relation between $E$, and $l$ and $e$

\begin{align*}
(Lv - k \hat{r})^2 &= k^2 e ^2 \\
(Lv)^2 - 2k(Lv) \cdot \hat{r} + k^2 &= k^2 e^2 \\
l^2 v^2 -2k \frac{\left< Lvr \right>}{r} &= k^2 (e^2 -1)
\end{align*}

Using $\eqref{eq3}$

\begin{align*}
l^2 v^2 - \frac{2kl^2}{mr} &= k^2 (e^2 -1)\\
\frac{2 l^2}{m} \left( \frac12  m v^2 - \frac kr \right) &= k^2 (e^2-1) \\
E &= \frac{m k^2}{2 l^2} (e^2 -1)
\end{align*}

Thursday, April 24, 2014

Kalman Filter

The Kalman filter describes the evolution of a state estimate with time. We're concerned with a state described by a probability distribution. The distribution is assumed to be Gaussian, so it's described with a mean $x$ and a variance $\Sigma$.

The system evolves according to a model $F$. So $x_1 = F x_0$ and $\Sigma_1 = F \Sigma_0 F^T$.

Additionally, we have a measurement of the state. This is also assumed to be a Gaussian distribution, with mean $z$ and variance $R$.

The Kalman trick is to realize that we have, at this point, two distributions describing our state knowledge from different avenues, and we should combine these to get our actual current knowledge distribution $x$ and $\Sigma$. For Gaussian distributions this is particularly simple, and the combined distribution is still Gaussian, so can be described with a combined mean and covariance.

$$x = x_1 + \frac{\Sigma_1 (z - x_1)}{\Sigma_1 + R}$$
$$\Sigma  = \Sigma_1 - \frac{\Sigma_1^2}{\Sigma_1 + R}$$

There's a common factor here, call it $K$

$$K \equiv \frac{\Sigma_1}{\Sigma_1 + R}$$

This is sometimes called the Kalman gain.

\begin{align}
x &= x_1 + K (z - x_1) \\
\Sigma &= \Sigma_1 - K \Sigma_1
\end{align}

From here, let us assume a multi-variate situation, so that $x$ is a vector, and use $\Sigma$ for the covariance matrix, instead of $\Sigma$ for variance.

Just one more thing gets us to a description of the Kalman update. The measurement space, with its particular units, etc., might not be the same as our state space. Let's use $H$ to get back and forth between these domains: $\hat{x} = H x $ and $\hat{\Sigma} = H \Sigma H^T$ We can now write the Kalman update:

\begin{align}
\hat{x}  &= \hat{x_1} + K (z-\hat{x_1}) \\
\hat{\Sigma} &= \hat{\Sigma_1} -K \hat{\Sigma_1} \\
K &\equiv  \hat{\Sigma_1} [ \hat{\Sigma_1} + R ]^{-1}
\end{align}

\begin{align}
H x  &= H x_1 + K (z-H x_1) \\
H \Sigma H^T &= H \Sigma_1 H^T -K H \Sigma_1 H^T \\
K &\equiv  H \Sigma_1 H^T [ H \Sigma_1 H^T  + R ]^{-1}
\end{align}

We can knock an $H$ off the front of the equation for $x$ and off the front and back of the equation for $\Sigma$, giving the final update relations.

\begin{align}
x  &= x_1 + K' (z-H x_1) \\
\Sigma  &= \Sigma_1  -K' H \Sigma_1 \\
K' &\equiv \Sigma_1 H^T [ H \Sigma_1 H^T  + R ]^{-1}
\end{align}

Thursday, November 14, 2013

GVD, Chirp, and Pulse Spreading

The Fourier transform of a gaussian pulse shape is also gaussian.  If $f(x) =  \exp(- (\frac{x}{a})^2)$, then we can develop the transform as

\begin{align}F(k) &= \frac{1}{2 \pi} \int f(x) e^{-i k x} {dx} \\
&= \frac{1}{2 \pi} \int e^{- ( \frac{x}{a})^2}   e^{-i k x} {dx} \\
&= \frac{1}{2 \pi} \int \exp \left[ - \left( \frac{x}{a} \right)^2 -i k x \right] {dx} \\
&=  \frac{a \sqrt{\pi}}{2 \pi} \exp \left( - \frac{a^2 k^2}{4} \right)\\
&=  \frac{a}{2 \sqrt{\pi}} \exp \left[ - \left(\frac{a k}{2} \right)^2 \right]
\end{align}

Consider now the time-evolution of this pulse, assuming this shape (in $x$), at $t=0$ and allowing each $k$ component to evolve according to $e^{-i w(k) t}$.  That is, evolve through linear media with $\omega \equiv \omega(k)$

$$
f(x,t) = \int  \frac{a}{2 \sqrt \pi} \exp \left[ - \left(\frac{a k}{2} \right)^2 \right] \exp (i k x - i \omega(k) t){dk}
$$

Let $\omega(k)  =  \omega_0 +  \omega_1 k  + \omega_2 k^2 + $ H.O.T.

\begin{align}
f(x,t) &= \frac{a}{2 \sqrt \pi}  \int  \exp \left( - \frac{a^2}{4} k^2  + i k x  - i \omega_0 t -i \omega_1 k t - i \omega_2 k^2 t \right)  {dk}\\

&= \frac{a e^{-i \omega_0 t}}{2 \sqrt \pi}  \int  \exp \left[ -  \left(\frac{a^2}{4} + i \omega_2 t \right)  k^2  + i (x   -  \omega_1  t) k  \right]  {dk} \\

&= \frac{a e^{-i \omega_0 t}}{2 \sqrt \pi} \sqrt \frac{ \pi}{ \left(\frac{a^2}{4} + i \omega_2 t \right)}   \exp \left( - \frac{( x   -  \omega_1  t) ^2}{4 \left(\frac{a^2}{4} + i \omega_2 t \right) }   \right) \\

&= \frac {a e^{-i \omega_0 t}}{\sqrt {{a^2} + 4i \omega_2 t }}   \exp \left( - \frac{( x   -  \omega_1
t)^2}{(a^2 + 4 i \omega_2 t ) } \right) \\

&= \frac { e^{-i \omega_0 t}}{\sqrt {1 + \frac{4i \omega_2 t} {a^2} }}   \exp \left( - \frac{( x   -  \omega_1  t)^2}{a^2 \left(1+ \frac{4 i \omega_2 t}{a^2} \right) } \right)

\end{align}

We can see the pulse moves with velocity $\omega_1$, the group velocity, and spreads with time according to a factor controlled by $\omega_2$, the group velocity dispersion.  To see this more clearly, let's look at the spectrum of a gaussian with a complex width parameter.

Since
$$
\int \exp(- ax^2) \exp (-bx) dx = \sqrt \frac{\pi}{a} \exp \left( \frac {b^2}{4 a }\right)
$$

We can write
$$
 \int \exp(- \Gamma k^2) \exp (ikx) dk = \sqrt \frac{\pi}{\Gamma} \exp \left( -\frac {x^2}{4 \Gamma } \right)
$$

Or
$$
\exp(- \Gamma k^2)   \iff   \sqrt \frac{\pi}{\Gamma} \exp \left( -\frac {x^2}{4 \Gamma } \right)
$$

If $\Gamma = a + ib$ then

\begin{align}
\exp(- (a+ib) k^2) &= \exp(- a k^2) \bbox[ border:2px solid yellow ] {\exp(-ib k^2)} \\
& \iff   \sqrt \frac{\pi}{(a+ib)} \exp \left( -\frac {x^2}{4 (a+ib) } \right) \\
&= \sqrt \frac{\pi}{(a+ib)} \bbox[ border:2px solid cyan ] {\exp \left( -\frac {a x^2}{4 (a^2+b^2) } \right)} \exp \left( i \frac {b x^2}{4 (a^2+b^2) } \right)
\end{align}
and we can see that the imaginary part of $\Gamma$ both controls the $\bbox[ border:2px solid yellow ] {\text{chirp}}$ in the $k$ domain, and influences $\bbox[ border:2px solid cyan ] {\text{width}}$ in the $x$ domain.

Tuesday, February 14, 2012

Partitions and Permutations

How many ways can you arrange 5 books on a shelf? Easy: $ 5 \times 4 \times 3 \times 2 \times 1 \equiv 5! $ This is a permutation.

Now, how many ways can you pull 4 books from a stack of 10? This is the same as partitioning the ten books into a stack on the shelf, and a stack in your hands.   An ambiguity arises: does the order of the books drawn, matter?  If no matter, then we want the number of ways all the books can be permuted, divided by the number of ways the books on the shelf can be permuted, divided again by the number of ways the books in hand can be permuted:

$$
\frac{10!}{4! \, 6!}
$$

This is sometimes referred to as a combination: the order of selection does not matter.  If the order of books drawn does matter,  then those permutations must not be divided out:  $ \frac {10!} {6!}$, and this kind of selection is called a permutation.

To generalize: the number of ways to permute N things is $N!$ and if some number of those things are identical, divide out those permutations:  $ \frac{N!}{ n_1! \, n_2! \, n_3! \, ... }$

What if we allow replacement?  Like with a combination lock, where we select 3 numbers but the same number can appear all three times.  This situation can not be understood as a partition, but this case is simple: if we have K slots that can take N different values, the number of possibilities is $N^K$.

If, however, we allow replacement but order does not matter, the situation becomes more difficult to count.  For that case, the trick is to observe that choosing, say, four from 10 things is the same as tossing four (indentical) things into 10 bins.  In other words, how many ways can you permute 4 X's with 9 partitions?   X|  |  |  |X|  |  |XX|  |  |




Thursday, August 18, 2011

Transform of a Product


Define $F$ such that
$$
f(x) = \frac{1}{\sqrt {2 \pi}}  \int{F(k) e^{ikx} {dk}}
$$

and therefore

$$
 F(k) = \frac{1}{\sqrt {2 \pi}} \int {f(x) e^{-ikx}}
$$

Then we have

$$
\int{f(x) g(x) {dx} }  =  \int{F(k)G(-k) {dk} }
$$

Since $ G^\star(-k) \iff g^*(x)$, we have

$$
\int{f(x) g^\star(x) {dx} } = \int {F(k) G^*(k)  {dk}}
$$


Since $g(x) e^{-i \alpha x}  \iff G(k + \alpha)$, a similar manipulation gives

\begin{align}
\int{ f(x) [g(x) e^{-ikx}] {dx} } &= \int{ F(k') G( - k'+ k) ) {dk'}} \\
&= \int{ F(k') G( k - k') ) {dk'}}
\end{align}

Sunday, July 10, 2011

Fourier Integral and Periodic Functions

Consider a function $f(t)$, sufficiently well-behaved such that

\begin{align}
f(t) = &\int_{-\infty}^{\infty} {F(\omega) e^{i \omega t} d{\omega }}\\
F(\omega) = \frac 1 {2 \pi} &\int_{-\infty}^{\infty}{f(t) e^{-i \omega t} dt} 
\end{align}

Periodic Function


What does the transform $F(\omega)$ look like for a periodic function?  Let $f(t)$ be periodic: $f(t)$ = $ \sum_n {f_0(t-nT) } $  where $f_0$ is a time-limited function defined over the interval $[0, T]$.

Compute $F(\omega)$:

\begin{align}

F(\omega) &= \frac 1 {2 \pi} \int_{-\infty}^{\infty}{ \sum_{n=-\infty}^{\infty} f_0(t-nT) e^{-i \omega t} {dt} }\\

 &= \sum_{n=-\infty}^{\infty} \frac 1 {2 \pi} \int_{-\infty}^{\infty}{  f_0(t-nT) e^{-i \omega t} {dt} } \\

 &= \sum_{n=-\infty}^{\infty} \frac 1 {2 \pi} \int_{-\infty}^{\infty}{  f_0(t) e^{-i \omega t}e^{-i \omega nT} {dt}} \\

 &= \bbox[ border:2px solid yellow ]{\left( \sum_{n=-\infty}^{\infty} e^{-i \omega nT} \right)}  \bbox[ border:2px solid cyan]{\left(   \frac 1 {2 \pi} \int_0^T{  f_0(t) e^{-i \omega t}{dt} } \right)}

\end{align}


The yellow box is  ${\Omega} \sum_{n=-\infty}^{\infty} \delta (\omega-n \Omega)$, where $\Omega = \frac {2 \pi} T$

So we end up with
\begin{align*}
F(\omega) &= \Omega \sum_{n=-\infty}^{\infty} \delta (\omega-n \Omega)  F_0(\omega) \\
&= \Omega \sum_{n=-\infty}^{\infty} \delta (\omega-n \Omega)  F_0(n \Omega)
\end{align*}

$F(\omega)$ is just a sampled version of the transform of $f_0(t)$.

Sampled Function


Working the other way now, let's transform a sampled function.

Let $F$ be the transform of $f$:

$$
F(k) =   \frac 1 {2 \pi} \int_{-\infty}^{\infty}{f(x) e^{-i k x} dx}
$$
$f(x)$ is a sampled version of $g(x)$: 
$$
f(x) = \sum_{n=-\infty}^{\infty} \delta(x - nL) g(x)
$$

The result is immediate:
$$
F(k) =  \frac 1 {2 \pi} \int_{-\infty}^{\infty}{  \left(  \sum_{n=-\infty}^{\infty} \delta(x - nL) \right) g(x)  e^{-i k x} dx}
$$

This looks like a transform of a product, which resolves to a convolution product of the individual transforms. Since $\sum_{n=-\infty}^{\infty} \delta(x - nL) \iff   \bbox[ border:2px solid yellow ]{\frac {\Omega} {2 \pi}  \sum_ {n=-\infty}^{\infty}\delta(k - n\Omega)}$, we can write

\begin{align}
F(k)  &=  \int_{-\infty}^{\infty} {  \left(  \frac {\Omega} {2 \pi} \sum_{n=-\infty}^{\infty} \delta(k - k^{'} - n\Omega) \right) G(k^{'})  dk^{'} } \label{ISI} \\
&=  \frac {\Omega} {2 \pi} \sum_{n=-\infty}^{\infty}  G(k - n \Omega)
\end{align}

This is a set of shifted copies of the transform of $g$, which might be a mess unless $g$ is zero outside the interval controlled by the sample rate $\Omega$.  Although, we will return to this later with an interesting example of when the overlap works out nicely.

Periodic, Sampled Function


Finally, let's work out the transform of a periodic, sampled function.

\begin{align}
F(\omega) &= \frac {1}{2 \pi} \int f(t) e^{-i \omega t} {dt} \\
&=\frac {1}{2 \pi} \int \sum_{m=-\infty}^{\infty} \sum_{n=0}^{L-1} d_n \delta(t-mP-nT) e^{-i \omega t} {dt}\\
&=\frac {1}{2 \pi}\sum_m \sum_{n=0}^{L-1} d_n e^{-i \omega m P} e^{-i \omega n T}\\
&=\frac {1}{2 \pi}\bbox[ border:2px solid yellow ]{\sum_m e^{-i \omega m P}} \sum_{n=0}^{L-1} d_n e^{-i \omega n T} \\
&= \frac {1}{P}\sum_m \delta(\omega - \frac {m 2 \pi}{P}) \sum_{n=0}^{L-1} d_n e^{-i \omega n T} \\
&= \sum_m \delta(\omega - \frac {m 2 \pi}{P}) \frac {1}{P} \sum_{n=0}^{L-1} d_n e^{-i 2 \pi m n  \frac {T}{P} } \\
&= \sum_m \delta(\omega -  m \Omega) \frac {\Omega}{2 \pi} \sum_{n=0}^{L-1} d_n e^{-i 2 \pi m \frac {n} L }  \\
&= \sum_m \delta(\omega -  m \Omega)  c_m
\end{align}

We used $\Omega P = 2 \pi$ and $P \equiv T L$, and defined

\begin{align*}
c_m &= \frac {\Omega}{2 \pi} \sum_{n=0}^{L-1} d_n e^{-i 2 \pi m \frac {n} L } \\
&= \left({1 \over T}\right)  \left({1 \over L} \sum_{n=0}^{L-1} d_n e^{-i 2 \pi m \frac {n} L }\right)
\end{align*}

In general, the transform of a periodic function is discrete, the transform of a discrete function is periodic, and the transform of a periodic and discrete function is periodic and discrete.

The Nysquist Criterion for Avoiding ISI


Returning to equation \eqref{ISI}, consider the following system

$$ y(t) = \sum_n c_n h(t-nT)   $$

Sampled at intervals $t-kT=0$, this looks like

\begin{align*}
y(kT) &= \sum_n c_n h(kT-nT) \\
\hat y(k) &= \sum_n c_n  \hat h(k-n)
\end{align*}

We would like to investigate the situation where $\hat y(k)$ is just $c_k$.  That is, $\hat h(k)$ are all zero with one exception: $\hat h(k) = \delta_k$.  So let's look at the Fourier analysis of a discrete version of $h(t)$:

\begin{align*}

\sum_{n=-\infty}^{\infty} \delta(t - nT) h(t) &= \int \frac {\Omega} {2 \pi} \sum_{n=-\infty}^{\infty}  H(\omega - n \Omega) e^{i \omega t} {d\omega} \\

\delta(t) &=  \int \frac {\Omega} {2 \pi} \sum_{n=-\infty}^{\infty}  H(\omega - n \Omega) e^{i \omega t} {d\omega}

\end{align*}

So we must have
$$
\Omega \sum_{n=-\infty}^{\infty}  H(\omega - n \Omega) = 1
$$

In the study of communication in bandlimited channels, this is known as the Nyquist Criterion for the elimination of inter-symbol interference.

Monday, May 2, 2011

Homomorphism

Group

Let $G$ be a group with elements ${g}$, and $H$ is a homomorphism of $G$. The set of elements mapped by $H$ to the identity is the kernel $K$ of $H$. $K$ is a normal subgroup of G, and the cosets of $K$ form another group ${G \over K}$, called a factor group.

Ring

Let $R$ be a ring, and $H$ is a homomorphism of $R$. The kernel $K$ of $H$ is an ideal of $R$, somewhat analogous to the normal subgroup discussed above. The cosets of the ideal $K$ form a factor ring (sometimes: quotient ring) ${R \over K}$. We would like to know how the properties of an ideal $K$ determine the properties of the associated factor ring. We will let our rings be commutative and have identity.

Elements that multiply to zero are called zero divisors. In $Z_6$ (the integers modulo 6), 2 and 3 multiply to zero, so they are zero divisors. If $n$ is not prime, $Z_n$ will have zero-divisors. A ring with no zero divisors is an integral domain. In an integral domain, $ab=0 \implies a=0$ or $b=0$

The elements that have inverses are called units.  The units form a group.  The group of units in $Z_6$ are {1,5}.

A ring in which every element has a multiplicative inverse is a field. A field is necessarily an integral domain. If $ab =0$ , and $a \neq 0$, then multiply by $a^{-1}$ and conclude $b=0$. That the integers $Z$ is not a field shows the reverse is not true. It is true, however, that every finite integral domain is a field.

Principal Ideal

The ideal $(b)$ consisting of all multiples of an element $b$ of $R$ is a principal ideal. This is the smallest ideal containing $b$. If it's unclear how an ideal might *not* be principal, consider a polynomial in two variables over the complex numbers. The ideal generated by x and y is not principal, because if it's principal there's a generator, and if there's a generator, say $p$, this divides every element of the ideal, but that must be a constant (non-zero). Yet there are no constants in the ideal, contradiction.

We can have an ideal generated by a subset, but this is a principal ideal only if the subset is a single element.

Prime Ideal

An ideal $K$ is prime if and only if $ab$ in $K$ implies $a$ or $b$ is in $K$. If $K$ is a prime ideal, the factor ring ${R \over K}$ is an integral domain.

Maximal Ideal

An ideal $K$ is maximal if the only ideal containing $K$ is $R$. If $K$ is a maximal ideal, ${R \over K}$ is a field. A maximal ideal must be prime, but the reverse is not true: a prime ideal might not be maximal.

Friday, April 16, 2010

Bias and Variance

Consider $\langle (y-\theta)^2\rangle $, the expectation of the squared error between estimator $y$, a random variable, and target $\theta$ (determined, not random).

After a small manipulation $ \langle (y-\langle y \rangle +\langle y\rangle-\theta)^2\rangle $, we can write

\[\langle(y-\langle y\rangle)^2+(\langle y\rangle-\theta)^2+2(y-\langle y\rangle)(\langle y\rangle-\theta)\rangle\]

\begin{equation*}
\langle (y-\langle y\rangle)^2\rangle+\langle (\langle y\rangle -\theta)^2\rangle+2 \langle y-\langle y\rangle \rangle \langle \langle y\rangle-\theta\rangle
\end{equation*}


But since $\langle y- \langle y\rangle \rangle$ is $0$, we can write

\begin{equation} \label{bv1}
\langle (y-\langle y\rangle)^2\rangle+ (\langle y\rangle-\theta)^2
\end{equation}

That is, the variance plus the bias squared. An immediate consequence is that if your estimator is unbiased, the mean squared error of your estimator is the same as your estimator's variance. That's a nice result.

Let's try to generalize a little. What if $\theta$ is a random variable? Say, $\theta = \hat{f} + \epsilon$, where $\epsilon$ is a gaussian random variable, center zero and width $\sigma^2$.

Beginning with \eqref{bv1}, we have

\[ \langle (y-\langle y\rangle)^2\rangle+\langle (\langle y\rangle -(\hat{f} + \epsilon))^2\rangle\]

When expanding the second term, all cross terms with $\epsilon$ vanish (center is zero), so we have

\[ \langle (y-\langle y\rangle)^2\rangle+(\langle y\rangle - \hat{f})^2 + \langle\epsilon^2\rangle \]

Or the variance plus the squared bias plus an unavoidable squared error term.

Friday, March 6, 2009

Lagrangian Description of the Electron-Photon Interaction

$\newcommand{\abs}[1] {\left\vert #1 \right\vert}$

The Lagrangian is described in three parts: $L = L_P + L_R + L_I $ (particles + radiation + interaction).  We use a Lagrangian density, so the action is $ S=\int L dt = \mathscr{\iint L \;dx\; dt} $

\begin{align*}
L_P &= \sum_\alpha {1 \over 2} m_\alpha \dot{x}_\alpha^2 \\
L_R &=  {\epsilon_0 \over 2}  \int{E^2 -c^2 B^2 dr}
\end{align*}


\begin{align*}
L_I &= \sum_\alpha  q_\alpha \dot{r}_\alpha \cdot A(r_\alpha) - q_\alpha U(r_\alpha) \\
&= \int  J(x) \cdot A(x) - \rho U(x) dx
\end{align*}

Rendering the integrals in Fourier space:
$\require{enclose}$

\begin{align*}
L_R &= \epsilon \, \enclose{horizontalstrike}{\int} \mathscr{ \abs{E}^2  - c^2 \abs{B}^2 } dk\\
L_I &= \enclose{horizontalstrike}\int  \mathscr{j^* \cdot  A + j \cdot A^*- \rho^* U- \rho U^*} dk
\end{align*}

But since $\mathscr{E =-i k U- \dot A}$ and $\mathscr{ B = i k \times A}$ we can write
$$ \mathscr{ \abs{E}^2 = k^2 \abs{U}^2 +  \abs{\dot A}^2 + i  U k \cdot \dot A^*} + \text{c.c.}$$

Notice that $\dot U$ does not appear in $L$

$$ 0 = \mathscr{ \frac{\partial L}{\partial U^*} = \epsilon_0 (k^2 U - i k \cdot \dot{A}) - \rho} $$
$$ \mathscr{ U = \frac{\rho}{\epsilon_0 k^2} + {i \over k} \dot{A_{||}}} $$

Use this to eliminate U:

\begin{equation*}
\mathscr{E = -{ i \rho \underline{k} \over \epsilon_0 k^2} + \dot {A_{||}} - \dot{A} = -{i \rho \underline{k} \over \epsilon_0 k^2} - \dot{A_{\perp}}}
\end{equation*}

Which lets us write
\begin{align*}
L_R &= \mathscr{ \epsilon_0 \enclose{horizontalstrike}\int {\abs{\rho}^2 \over k^2 \epsilon_0^2} + \abs{\dot{A_{\perp}}}^2 -c^2 k^2 \abs{A_{\perp}}^2 dk} \\
L_I &= \mathscr{ \enclose{horizontalstrike}\int j_{||} A_{||}^* + \text{c.c.} - {2 \abs{\rho}^2 \over \epsilon_0 k^2} -\frac{i \rho^* \dot{A_{||}}}{k} + \text{c.c.} \;  dk}
\end{align*}

Observe that $ \mathscr{A}_{||}$ only appears in $\mathscr{L}_I$

\begin{align*}
\mathscr{\frac{\partial L}  {\partial A_{||}^*} }  &= j_{||}     \\
\mathscr{\frac{\partial L}{ \partial \dot{A_{||}^*}}} &= \frac{i \rho}{k}  \\
j_{||} &= {i \dot{\rho} \over k}
\end{align*}

Use this to eliminate $j_{||}$

\begin{align*}
L_I &= \mathscr{  \enclose{horizontalstrike}\int i {\dot{\rho} \over k} A_{||}^* + \text{c.c.}  - 2 {\abs{\rho}^2 \over \epsilon_0 k^2} - i {\rho^* \dot{A}_{||} \over k} + \text{c.c.} +  j_{\perp} \cdot A_{\perp}^* + \text{c.c.} \;  dk  }\\
&= \mathscr{  \enclose{horizontalstrike}\int  -2 {\abs{\rho}^2 \over \epsilon_0 k^2}  + {i \over k} {d \over dt}  (\rho A_{||}^* - \text{c.c.}) +  j_{\perp} \cdot A_{\perp}^* + \text{c.c.} \; dk }
\end{align*}

\begin{align*}
L = \sum_\alpha \frac12 m_\alpha \abs {\dot{r}_{\alpha}}^2 &- \enclose{horizontalstrike}\int {\abs{\rho}^2 \over \epsilon_0 k^2}  dk \\
&+ \mathscr{ \epsilon_0 \enclose{horizontalstrike}\int \abs{\dot{A_{\perp}}}^2 - c^2 k^2 \abs{A_{\perp}}^2 dk }\\
&+ \mathscr{ \enclose{horizontalstrike} \int j_{\perp} \cdot A_{\perp}^* + \text{c.c.} \;  dk}\\
&+ \mathscr{ {i \over k} {d \over dt} \enclose{horizontalstrike}\int \rho A_{||}^* -\text{c.c.} \;dk }
\end{align*}

Lagrangian dynamics is insensitive to a total time derivative, so the last term gives some freedom in choice of $A_{||}$ ("gauge").


Monday, July 7, 2008

Birthday Problem

You have N people at a dinner party.  How likely is it that two guests have the same birthday?

The first insight is that the probability $A$ of having a birthday pair is one, less the probability $B$ of having $N$ totally unique birthdays. That is, $A=1-B$.  Let's work $B$, it makes for easy counting.

We're looking for a probability, which in this case is a ratio between two things.  The first thing we want to know is, how many ways are there for $N$ people to have birthdays (any birthdays at all)?  That's easy: $365^N$.  The other thing we want to know is, how many ways are there for N people to NOT share a birthday?  That is, how many ways can we have $N$ unique birthdays?  Well, for more than 365 people we must have at least one pair, but for less than this a little reflection reveals the answer to be $365 * 364 * 363 * 362 * ... * (365-N)$.  In other words, let one person have any birthday, the next can have almost any birthday, the next one less, etc... The ratio of these is the probability we're after.

$$ P(N) = \frac{365 * 364 * 363 * ... * (365-N)} {365^N} =  \frac{365!}{365^N (365-N)!} $$

Here's a plot of P(N) for up to 100 people.  You can see the probability of unique birthdays is very high for small numbers of people, but crosses the "even odds" plane at 23 people.

Friday, June 13, 2008

Amort

Ever thought about what the actual formula is for your loan payment schedule? Spreadsheet functions and financial calculators make it easy to avoid the details, but let's dig into it anyway. The formula

\[ A = P \left( \frac{r  }{1-(1+r)^{-N}} \right) \]

looks sorta complicated.

Let $A$ be the monthly payment we're after, and $P$ is the initial principle. The monthly interest rate is $r$.  That is, the loan grows by a factor of $(1+r)$ every month. We want to set up the payments so that after N months, the loan is paid and we're done.

So how do we do this? Every month the balance on the loan will be different, let's use $\alpha_m$ for the balance at month $m$. We know two things right away:

\[ \alpha_{0} = P \tag{initial balance} \]
\[ \alpha_{N} = 0 \tag{ending balance}\]

So we have

\[ 0 = \alpha_N  = (1+r) \alpha_{N-1} - A \]

The last balance, $0$, is the penultimate balance increased by the interest and reduced by the monthly payment $A$. Let's rewrite in terms of $A$:

\[ A = (1+r) (\alpha_{N-1}) \]

\[ A = (1+r)( (1+r) \alpha_{N-2} - A ) \]

 \[ A (1  + (1+r))  = (1+r)^2 \alpha_{N-2} \]

 \[ A (1  + (1+r))  = (1+r)^2 ((1+r) \alpha_{N-3}  -A)\]
 
  \[ A (1  + (1+r) + (1+r)^2)  = (1+r)^3 \alpha_{N-3} \]

We can take this all the way to N, where the form for $A$ will be:

\[ A \left(  \sum_{n=0}^{N-1}  (1+r)^n \right) = (1+r)^N \alpha_0 \]

 The sum on the left is a geometric series that we know how to evaluate:

\[ A \left(  \frac{1- (1+r)^N}{1-(1+r)} \right) = (1+r)^N P  \]

Finally,

\[ A = \frac{rP(1+r)^N } {(1+r)^N -1}= P \left( \frac{r} {1-(1+r)^{-N}} \right)  \]

Saturday, June 9, 2007

Divergent Series

Consider these sums:

\begin{equation} \label{eq1}
\sum_{n=0}^\infty {(1)^n} = 1+1+1+1+1+ \cdots
\end{equation}

\begin{equation} \label{eq2}
\sum_{n=0}^\infty {(-1)^n} = 1-1+1-1+1-1+ \cdots
\end{equation}

\begin{equation} \label{eq3}
\sum_{n=0}^\infty {2^n} = 1+2+4+8+16+ \cdots
\end{equation}

Is there a way to give a meaning to these sums?  That is, to sum them?

What if there was?  Imagine that we have a machine $ \mathfrak{S}$ that can sum infinite numbers of numbers. What would such a thing give for these sums? Let's consider what basic sanity we might expect. To start, it seems like anything that adds should maintain this kind of consistency:

\begin{equation*}
\mathfrak{S} (a+b+c+d+ \cdots) = a+ \mathfrak{S} (b+c+d + \cdots) \tag{property 1}
\end{equation*}

This is sometimes called stability. In addition, let's require our machine can maintain linear behavior:

\begin{equation*}
\mathfrak{S} \left[ \sum_n ( \alpha  a_n + \beta  b_n) \right] = \alpha \, \mathfrak{S} \left( \sum_n a_n \right) + \beta \, \mathfrak{S} \left( \sum_n b_n \right)   \tag{property 2}
\end{equation*}

Let's apply our machine to \eqref{eq1}: $ s  = \mathfrak{S} (1+1+1+ \cdots) $

Using property 1, we can write $ s  = 1 + \mathfrak{S} (1+1+1+ \cdots) $. So then $ s = 1 + s $. Seems ridiculous, unless $ s $ is infinity, which is what we might expect for this sum. So maybe things are working out so far.

Let's look at \eqref{eq2}.

\[ s = \mathfrak{S}(1-1+1-1+1-1+\cdots) \]

\[ s = 1+ \mathfrak{S}(-1+1-1+1-1+\cdots) \]

\[ s = 1- \mathfrak{S}(1-1+1-1+1-1+\cdots) \]

So then $ s = 1-s $, or $ s=1/2 $.  Well that's disappointing.  The individual terms of the sum bounce around interminably, the partial sums do the same, but our $ \mathfrak{S}$ gives a small positive rational number.  This is getting weird.  Did we go wrong somewhere?  The consistency properties we assumed were pretty modest.

Let's look at \eqref{eq3}:
\[ \sum_{n=0}^\infty 2^n = 1 +2 +4  +8 +16 + \cdots \]

\[ s = 1 + \mathfrak{S}(2+4+8+16+ \cdots) \]

\[ s = 1 + 2  \mathfrak{S}(1+2+4+8+ \cdots) \]

So $ s= 1 + 2s $, and $s =-1$.

Wow.  The terms of this sum are, each, bigger than in \eqref{eq2}, but now we are getting a finite number and it's negative.  Let's try one more.

\[ s = \mathfrak{S} (1 -2 +3 -4 +5 -6 +7 -8 + \cdots ) \]

\[ s= 1+ \mathfrak{S}(-2+3-4+5-6+7-8 + \cdots) \]

Add these two together (we're using property 2 this time):

\[ 2s = 1 + \mathfrak{S}(-1+1-1+1-1+1-1+ \cdots) \]

\[ 2s = 1 - \mathfrak{S}(1-1+1-1+1-1+1 -1 +\cdots) \]

From our result for \eqref{eq2}, we can write $ 2s = 1 - \frac{1}{2}$, or $ s = \frac{1}{4} $.

All of these sums are divergent -- the partial sums do not converge.  Perhaps divergent sums are just trouble and we're better off not trying to find meaning in that demonic mess.  Let's back up to safer ground and examine an infinite sum that behaves nicely and converges.

This sum

\begin{equation*}
\sum_{n=1}^\infty {\frac{1}{n}} \tag{harmonic}
\end{equation*}

(famously) diverges.  But this one

\begin{equation*}
\sum_{n=1}^\infty {\frac{(-1)^{n-1}}{n}} \tag{alternating harmonic}
\end{equation*}

converges, very slowly, to $ ln(2)$. After a moment's consideration you may see that the alternating harmonic series contains two sums, in a sense.  One has an infinite number of positive terms, $ 1 + \frac{1}{3} + \frac{1}{5} +\cdots$ and the other has an infinite number of negative terms $ - \frac{1}{2} - \frac{1}{4} - \frac{1}{6} -\cdots$. So if we are allowed to associate terms at will, I can make the sum converge to any number I want.  Say, 50.  I just sum terms from the positive set until I get just over 50, then start using terms from the negative set till I get under 50, then use positive terms...

So it should be clear that we can't insist on associativity when working with infinite sums.  That is, addition is not infinitely associative. This is surprising. If readers were wondering why we didn't allow the infinite summing machine $ \mathfrak{S}$ to be associative, now you know.

Okay then, it seems unavoidable that infinite sums are weird, and convergence doesn't relieve the issue.  Something spooky is going on.  Can we come up with a concrete example of something that behaves like $ \mathfrak{S} $?  The following approach, due to Euler, might illuminate things.

Euler Summation

Consider an infinite sum $ s = \sum_n a_n$, and form a function $ f(x) $ like so:

\[  f(x) = \sum_n a_n x^n \]

The Euler sum is

\[ \lim_{x \to 1-}f(x) \]

Let's try this on \eqref{eq2}.

\[ Euler(1-1+1-1+1-1+\cdots) = 1-x + x^2-x^3 + \cdots = \lim_{x \to 1-}\frac{1}{1+x} = \frac{1}{2} \]

Well that's nice, it's consistent with our previous approach. But maybe we're just lucky.  Let's try something else.

Borel Summation

Consider $ \int_0^\infty \frac{t^n}{e^t}dt = n!$, so that

\[ 1= \int_0^\infty \frac{t^n}{{e^t}{n!}}dt\]

Form a new series
\[ \sum_n a_n {\int_0^\infty \frac{t^n}{{e^t}{n!}}dt}\]

The Borel sum is what you get if you reverse the order of summation and integration:

\[ \int_0^\infty dt \, e^{-t} \sum_n \frac{t^n a_n}{n!}\]

Now look at $ Borel(1-1+1-1+1-1 + \cdots)$.

\[ \int_0^\infty dt \, e^{-t} \sum_n \frac{{(-t)}^n}{n!} \]

\[ \int_0^\infty e^{-2t} dt = \frac{1}{2}\]

Same!

General Summation

Returning to our general summation approach, let's investigate the impact of this slight modification to \eqref{eq2}:
\[ s = \mathfrak{S} (1 +0 -1 +1 +0 -1 +1 +0 -1 + \cdots) \]
\[ s = 1+ \mathfrak{S} (0 -1 +1 +0 -1 +1 +0 -1 + \cdots) \]
\[ s = 1+ \mathfrak{S} (-1 +1 +0 -1 +1 +0 -1 + \cdots) \]
Adding these three relations together: $ 3s=2 +  \mathfrak{S}(0+0+0+0+0+\cdots) $, or $ s=2/3 $.  So inserting zeros into the sum changes the result from what we got for \eqref{eq2}.  It's hard to rely on intuition in this area.  In fact, we might now question what really constitutes the "basic sanity" we started with.

Zeta Function Regularization

Consider the Zeta function:  $ \zeta(z) = \sum_{n=1}^\infty \frac{1}{n^z}$.  Zeta is remarkable for many reasons and the first is that it is a meromorphic function with one and only one simple pole, at $ z=1 $.  At that point, $\zeta(1)$ is the harmonic series we met earlier, which we were not able to assign a sum to. But look what $\zeta$ says about $z=0$:

\[\zeta(0) = \sum_{n=1}^\infty \frac{1}{1} = 1+1+1+1+1+1+1+ \cdots \]

$ \zeta(z)$ is analytic at $z=0$, and in fact $\zeta(0) = - \frac{1}{2} $.  What does this mean, especially after the happy result we got from general summation?  Well, $\zeta(z)$ must be either not stable or not linear, or both.  Whatever it is, that a perfectly analytic function assigns a value to this sum is not easy to dismiss.

There is a whole field of study that goes under the name summation theory that explores summation strategies and axioms.