Generalized coordinates

In this post I’ll try to explain what Generalized coordinates are, how they came to be and finally where they are used with a simple.

To start with imagine a system of $k$ point masses present in a $d$ dimensional space. The $i^{th}$ point has a position vector $r_i = (x_{i1}, x_{i2}, \cdots, x_{id})$.

Note that if the system was unconstrained, then as any of the $k$ bodies can occupy any point in the space, we will need $kd$ number of scalars to completely define any particular configuration of the system.

For visualization purposes consider $d=3$, so we use $3$ scalars (the $x, y, z$ coordinates) to define the position of each point mass. As we have a total of $k$ such masses we need $3k$ scalars.

But often the systems we need to analyse are under some sort of constraints. In case of a robotic arm, each hinge point always maintains a fixed distance between them due to the rigid rod connecting them.

In such situations we don’t need all $kd$ scalars, rather we can use less number of scalars and still manage to uniquely represent all possible states that this constrained system can take. This is because the constraints decrease the degrees of freedom available to that system.

Holonomic constraints

We classify constraints as either being holonomic or nonholonomic.

If a constraint can be expressed purely as a relationship between the position variables ($r_1, r_2, \cdots, r_k$) and time $t$.
$${\displaystyle f(r_1, r_2, \cdots, r_k,\ t)=0,\,}$$
Then using each such constraint equation we can reduce the number of scalars needed by one. If we knew the time $t$ and all but one position variable, using the constraint equation that missing coordinate can be deduced.

When the constraint has inequalities or higher order derivative terms in them, they are classified as nonholonomic. Note that we have an equality symbol in the above expression. Also the function has to of the position variables themselves, not their derivatives.

All constraints which are unable to be expressed in the above form are called nonholonomic constraints.

Note that if there were velocity terms ($\dot r_i$) in the constraint equation, then integrating that constraint could in some cases yield an equation of the form shown above, if so then it’s a holonomic constraint. So what Wikipedia means when they say that if the constraint has velocity terms they “are not usually holonomic” is that it’s kinda hard to bring it to this form in that case.

When constraints are purely dependent on the positions at any time $t$ its a holonomic constraint, but in a nonholonomic constraint the constraint changes depending on how you reach that position! This seems to be analogous to path and state functions.

To visualise how a holonomic constraint reduces DOF, imagine a point mass is constrained to move along the line $y = mx+c$. So without this constraint we needed both position variables $(x,y)$ to describe its position anywhere in the XY plane, but now we just need $(x)$. As if we have $x$, its position is now known to be $(x,mx+c)$ due to the constraint.

Constrained to move in a circle of radius $r$

So a point mass is moving around in a plane (say the $XY$ plane). We can describe its position at any point of time using the vector,

$$\vec r = \left [ \begin{array}{cc} r_1(t) \\ r_2(t) \\ \end{array} \right ]= \left [ \begin{array}{cc} x(t) \\ y(t) \\ \end{array} \right ]$$

If we are given that this point mass is constrained to move in a circle of radius $c$, i.e. the constraint equation would like this,

$$|\vec r| = c^2$$
$$x(t)^2 + y(t)^2 = c^2$$

Note how the equation can be written as $x(t)^2 + y(t)^2 - c^2 + 0t = 0$, i.e. in the form $f(x,y,t)=0$. Therefore we conclude that this is a holonomic constraint.

So the “generalized coordinate” I will use is $\theta$, the angle between the $x$ axis and the position vector $\vec r$.

We know that $x = c \cos \theta$ and $y = c \sin \theta$.

Let the point mass move a little, it can of course only move in a small arc along the circle due to the constraint imposed on it. Such displacement is called Virtual displacement. In terms of $\theta$ it’s $\delta \theta$, while in terms of $\vec r$ it’s

$$\delta\vec r = \left [ \begin{array}{cc} \delta x \\ \delta y \\ \end{array} \right ]$$

The following formula relates this change in the $m$ generalized coordinates ($\delta q$) to that of the change in natural coordinates ($\delta\vec r$).
$$\delta\vec r_i=\sum _{{j=1}}^{m}{\frac {\partial {\vec r}_{i}}{\partial q_{j}}}\delta q_{j}\,$$

The term ${\frac {\partial {\vec r}_{i}}{\partial q_{j}}}$ is a measure of how much change in ${\vec r}_{i}$ happens for a unit change in $q_{j}$, keeping all other generalized coordinates $q_i \ \forall \ i \in [1,... , j-1, j+1, ... m]$ constant - that’s what partial differentiation is.

So by multiplying ${\frac {\partial {\vec r}_{i}}{\partial q_{j}}}$ with $\delta q_{j}$, you get the contributed of the $j_{th}$ generalized coordinate - $\ q_j$ in the displacement $\delta\vec r_i$.

Summing over all these contributions (from $1$ to $m$) we get the net virtual displacement $\delta\vec r_i$.

Coming back to our example, we can therefore note that

$$\delta x = \frac{\partial x}{\partial \theta}\delta \theta = c(-\sin \theta)\delta \theta$$

$$\delta y = \frac{\partial y}{\partial \theta}\delta \theta = c(\cos \theta)\delta \theta$$

Credits to Maschen - for this image, CC0, Link

This makes sense as for small displacement, the $\delta\vec r$ shown in the figure will be a straight line, thus instead of a sector, we have a right angled triangle with $\delta x$, $\delta y$ and $|\delta\vec r|$ as the sides. $|\delta\vec r|$ will be the length of the arc which is $c \ \delta \theta$.

So from pythagoras theorem we have $(\delta x)^2 + (\delta y)^2 = |\delta\vec r|^2 = (c \ \delta \theta)^2$.

And sure enough this fits in with our “derived” relationship between the generalized and natural coordinates.

$$(\delta x)^2 + (\delta y)^2 = (c(-\sin \theta)\delta \theta)^2 + (c(\cos \theta)\delta \theta)^2$$

$$=c^2(\delta \theta)^2((\sin \theta)^2+(\cos \theta)^2) = (c \ \delta \theta)^2$$

You can also try doing this with $q = s$ rather than $q = \theta$ like how I showed above. Where $s$ is the length of the curve traced out by the point mass from some reference point.

I thought of writing it up but it turns out that the length of a curve $y=f(x)$ from $(a,f(a)$ to $(b,f(b))$ is $S$. Which is kinda complicated,

$$S = \int_{a}^b\sqrt {1 + (f'(x))^2}$$

Well as always let me know if you have any thoughts on my post.

Written with StackEdit.

Dimensionality Reduction

Dimensionality reduction is a method of reducing the number of random variables under consideration via obtaining a set “uncorrelated” principle variables.

If we have $n$ features (each of them are the random variables we just talked about) then each data sample will be a $n$ dimensional vector in the $\mathcal{R}^n$ feature space. when $n$ is a large number, patterns in data can be hard to find as graphical representation is not possible.

What if we are able to identify a set of “principal variables” which are transformations of the existing $n$ variables such that we don’t lose much information.

PCA is a famous example when such data transformation are linear, there are other non linear methods too.

Here I’ll be now talking about PCA (Principal Component Analysis) and K-LDA (kernelized version of linear discriminant analysis (LDA)] ) and using these techniques in order to reduce the 2 dimensional data in this dataset into a single dimensions.

Principal Component Analysis

PCA is the orthogonal projection of the data onto a lower dimensional linear space. such that the variance of the projected data is maximized.

Variance is spread of the data, if the data is more spread, the better it is for us as we can easily separate them into clusters etc. We lack information if all the data vectors are localized around the mean.

So to perform PCA, we calculate the directions (in the $n$ dimensional feature space) in which variance is maximum (say $m$ such directions). Using these directions as basis vectors we can now project all the data vectors onto this subspace.

Here $m<n$, and we choose the directions by calculating the eigenvectors of the covariance matrix. The eigenvectors corresponding to the larger eigenvalues will be the directions where the variance varies more.

Now I will be doing the above method in MATLAB, as part of assignment as I am in ICE department I am assigned the 6th dataset - TRAIN{1,6} so as you can see my - has 4 clusters -

• So we can load the data set and initialize the clusters. After which we can plot them as a scatter graph.

load TRAINTEST2D

cluster1 = TRAIN{1,6}{1,1}; % Green
cluster2 = TRAIN{1,6}{1,2}; % Blue
cluster3 = TRAIN{1,6}{1,3}; % Red
cluster4 = TRAIN{1,6}{1,4}; % Cyan


scatter(cluster1(1,:), cluster1(2,:), 'g'); 
hold on;
scatter(cluster2(1,:), cluster2(2,:), 'b');
hold on;
scatter(cluster3(1,:), cluster3(2,:), 'r');
hold on;
scatter(cluster4(1,:), cluster4(2,:), 'c');
hold on;

• So we have data comprising a set of observations of $p = 2$ variables. This data can be arranged as a set of $n = 52$ data vectors $\mathbf{x}_1 \ldots \mathbf{x}_n$ with each $\mathbf{x}_i$ representing a single grouped observation of the p variables.
  Here 4 clusters with 13 data vectors each, contributed to form $13 \times 4 = 52$. Write $\mathbf{x}_1 \ldots \mathbf{x}_n$ as row vectors, each of which has p columns. Place these row vectors into a single matrix X of dimensions $n \times p$.
• Now we center our data at the origin. This is supposed to help in reducing the Mean Square Error when we approximate the data. This answer also provides some insight into how it’s necessary to compute the covariance matrix. Note that that after centering at origin, if you execute sum(B) you will get a very small number close to 0. So we can conclude the data is successfully centered.
  But the images before and after centreing don’t look that different for this dataset as the empirical mean is $-0.1502 , -0.0465$ the entire data set just shifts as a whole fractionally.

% number of data vectors - 13*4
n = 52;
% number of dimensions (before reduction)
p = 2; 
% arrange data vectors as rows
X = [cluster1';cluster2'; cluster3'; cluster4'];
% find both x and y's mean - empirical mean vector
u = (1/n)*sum(X);
% subtract from the mean - Thus achieving centering
B = X - ones(52,1)*u;


% Plot the data after centering it around the Origin
figure(2);
scatter(B(1:13,1), B(1:13,2), 'g'); 
hold on;
scatter(B(14:26,1), B(14:26,2), 'b');
hold on;
scatter(B(27:39,1), B(27:39,2), 'r');
hold on;
scatter(B(40:52,1), B(40:52,2), 'c');
hold on;
legend('cluster 1','cluster 2','cluster 3','cluster 4');

• Next I will find the empirical covariance matrix from the outer product of matrix B with itself.

C = (B'*B)/(n-1);

• Then we can find the eigenvectors of the covariance matrix, choose the one which has a larger eigenvalue (In this data set $1.2>1.01$ ). Now we can project all the 2D points on the line (in the direction of that particular eigenvector).

% we get a column vector of eigen values - the first will be the one with
% the largest magnitude
E = svds(C);
% V will have the eigen vectors
[V,D] = eig(C);

% Thus the axis along which variance is Maximum is
maxVarianceAxis = V(:,2);

% No need to divide by (maxVarianceAxis'*maxVarianceAxis) as
% eig generates orthonormal eigenvectors
projectionMatrix = (maxVarianceAxis*maxVarianceAxis');

projCluster1 = projectionMatrix*cluster1;
projCluster2 = projectionMatrix*cluster2;
projCluster3 = projectionMatrix*cluster3;
projCluster4 = projectionMatrix*cluster4;

% Plot the data after projecting onto the principle component
figure(3);
scatter(projCluster1(1,:), projCluster1(2,:), 'g'); 
hold on;
scatter(projCluster2(1,:), projCluster2(2,:), 'b');
hold on;
scatter(projCluster3(1,:), projCluster3(2,:), 'r');
hold on;
scatter(projCluster4(1,:), projCluster4(2,:), 'c');
hold on;
legend('cluster 1','cluster 2','cluster 3','cluster 4');

As we can see the separation of blue and green clusters is good, but the light blue cluster is getting mixed up with the red cluster.

Before I wrap up, I’d like to show how the projection these data vectors onto the other eigenvector looks like,

See how the red and light blue is well separated here (But not as well separated green and the dark blue in the projection onto the eigenvector corresponding to the largest eigenvalue. The variance shown here is only slightly larger compared to the answer. (because the eigenvalues are close - about 0.2 difference) .

So here we had 2 features which we reduced to a single feature (the position of the data vector on the line - a single value if expressed as the distance from the origin). Thus we successfully implemented PCA and extracted a abstract feature where the data has maximum variance - globally thereby reducing dimensionality.

This leads us to LDA, where we no longer look at this same dataset as $52$ data vectors rather we take into consideration that that they are $13$ vectors from $4$ different classes.

Linear Discriminant Analysis

This is a method to find a linear combination of features that characterizes or separates two or more classes of objects.

This method is also known as Fisher’s linear discriminant, the idea is that we don’t wish to project such that the overall variance is maximized, rather we want the the clusters to be spread out (maximize the between-class covariance matrix after projection) while the data vectors of a particular class should be closer to the mean of that class (within-class covariance matrix for all the classes should be minimized).

The idea is that this enables us to easily discriminate new data vectors into classes. The more clearly separated the classes are, the less ambiguity we face when classifying.

I’ll explain using 2 classes and then extend it to multiple classes (in this dataset we are using 4 classes). let there be 2 classes - $C_1, C_2$, having $n_1, n_2$ number of vectors respectively. $x_{ij}$ is the $i^{th}$ vector in class $j$. So the mean of the classes will be,
$$m_1 = \frac{1}{n_1}\sum_{i = 1}^{n_1}x_{i1}$$
$$m_2 = \frac{1}{n_2}\sum_{i = 1}^{n_2}x_{i2}$$

the projection from $n$ to $1$ dimensions is done using $\vec w$,
$$y_{ij} = \vec w^Tx_{ij}$$
where the sizes are,
$$1 \times 1 = (1 \times n) \times (n \times 1)$$

The simplest measure of the separation of the classes, when projected along the line defined by $\vec w$, is the square of the separation between the projected class means. If the mean after projection of class 1 and 2 are - $\bar y_1$ and $\bar y_2$ then,
$$\bar y_2 - \bar y_1 = \vec w^Tm_2 - \vec w^Tm_1 = \vec w^T(m_2 - m_1)$$
$$(\bar y_2 - \bar y_1)^2 = \vec w^T(m_2 - m_1)(m_2 - m_1)^T\vec w = \vec w^TS_B\vec w$$

Here we define the between-class covariance matrix $S_B$ as,
$$S_B = (m_2 - m_1)(m_2 - m_1)^T$$

Thus $\vec w^TS_B\vec w$ is to be maximized.

While we minimize the sum of the within-class variance’s for each of the classes - $\sum_{j=1}^{k} \sigma^2_j$.

The within-class variance of the transformed data for the class $C_j$ is given by,
$$\sigma^2_j = \sum_{i = 1}^{n_j}(y_{ij} - \bar y_j)^2 = \sum_{i = 1}^{n_j}(\vec w^Tx_{ij} - \vec w^Tm_j)^2$$
$$\sigma^2_j = \sum_{i =1}^{n_j} \vec w^T(x_{ij} - m_j)(x_{ij} - m_j)^T \vec w$$
So after taking $\vec w$ terms out of the sigma and rewriting the summation in terms of the total within-class covariance matrix - $S_W$,

$$S_W = \sum_{i = 1}^{n_1} (x_{i1} - m_1)(x_{i1} - m_1)^T + \sum_{i = 1}^{n_2} (x_{i2} - m_2)(x_{i2} - m_2)^T$$

For 2 classes the measure of how close the data vectors are to thier respective classes mean is $\sigma^2_1 +\sigma^2_2$

$$\sigma^2_1 +\sigma^2_2 = \vec w^T S_W \vec w$$

Thus $\vec w^TS_W\vec w$ is to be minimized.

So with this we can construct the utility function which we need to maximize,

$$J(\vec w) = \frac{ |\vec w^TS_B \vec w|}{ |\vec w^TS_W\vec w|}$$

Differentiating $J(\vec w$) with respect to $\vec w$, setting it equal to zero, and rearranging gives,
$$\vec w \propto \mathbf{S}^{-1}_W(\mathbf{m}_2-\mathbf{m}_1)$$

Multi Class LDA

Thus we have found the direction of $\vec w$ for which our conditions are optimally achieved.

But extending this idea to multiple classes requires us to make some assumptions - so as to make the mathematics a little more easier.

Suppose that now there are $k$ classes. The within-class covariance matrix which is calculated is,
$$S_W = S_1 + ... + ... + S_n = \sum_{i = 1}^k \sum_{x\in c_i} (x-\bar x_i)(x-\bar x_i)^T$$

While the between-class covariance matrix is given by,
$$S_B = \sum_{i=1}^k n_i (\bar x_i - \bar x)(\bar x_i - \bar x)^T$$

In the above 2 equations $\bar x$ is the mean of all the data vectors, $\bar x_i$ is the mean of the data vectors in the $i^{th}$ class.

Using the revised covariance matrices and plugging them into the utility function. We will find that the optimal solution for $\vec w$ is,

$$S_W^{-1}S_B\vec w = \lambda \vec w$$

Optimal $\vec w$ is an eigenvector of $S_W^{-1} S_B$.

Kernel Trick + LDA = K-LDA

K-LDA the kernelized version of linear discriminant analysis (LDA) (also known as kernel Fisher discriminant analysis or kernel discriminant analysis) is when LDA is implicitly performed in a new feature space, which allows non-linear mappings to be learned.

For most real-world data a linear discriminant (offered by LDA) is not complex enough, because the data vectors might not be linearly separable at all.

We first map the data non-linearly into some feature space $\mathcal{F}$ (a higher dimensional space) and compute LDA there (For different kernel functions these mappings differ). The linear separation which LDA provides in the feature space will yield a non-linear discriminant in the original input space. You can think of a 2 clusters (in the 2d plane) with one of the clusters embedded inside the other, now if we map these vectors to 3 dimensions we can easily introduce a dummy variable (the third axis) such that the 2 clusters are easily separated linearly (by a plane). As you can see below the left side image is the Input space and the right side image is the feature space.

So as we are operating in the feature space, the only change in the LDA explained above is that we use $\phi (x_i)$ instead of $x$ everywhere. Explicitly computing the mappings $\phi(x_i)$ and then performing LDA can be computationally expensive, and in many cases intractable. For example, $\mathcal{F}$ may be infinitely dimensional (as in the case of the gaussian kernel).

So instead of designing a nonlinear transformation and calculating the transformed data vectors in the feature space - $\mathcal{F}$, we can rewrite the method in terms of dot products and using the kernel trick in which the dot product in the new feature space is replaced by a kernel function.

$$k(\mathbf{x},\mathbf{y})=\phi(\mathbf{x})\cdot\phi(\mathbf{y})$$

I used this paper and also refered from wikipedia - we can rewrite the equation for $J(\vec w)$ in terms of only dot products,

In the following I am assuming these sizes apply,
$$(\mathbf{x}_i^j)_{n \times 1} , (\alpha_i)_{1 \times 1} , (\vec w)_{n \times 1} , (\alpha)_{n_T \times 1} (m_i^{\phi} )_{n \times 1} , (M_i)_{n_T \times 1}$$
$$((M_i)_j)_{1 \times 1} , (\mathbf{K}_i)_{n_T \times n_i} , M_{n_T \times n_T} , N_{n_T \times n_T}$$

The remaining can be understood implicitly, here there are $n$ features, and the projection is onto 1 dimensions. These $n_T$ data vectors are separated into $k$ classes with the $i^{th}$ class having $n_i$ number of data vectors.

First of all $\vec w$ should lie on the span of all training samples $\mathcal{F}$ (read on theory of reproducing kernels to understand why)
$$\vec w = \sum_{i=1}^{n_T}\alpha_i\phi(\mathbf{x}_i)$$

where $n_T = \sum n_i$, is the total number of vectors.

Now the dot product of the mean of the $i^{th}$ class (in the feature space) $m_i^{\phi}$ with $\vec w$ will give,

$$\vec w^Tm_i^{\phi} = \alpha^T M_i$$

Here $\alpha$ is the set of $\alpha_i$ all piled together to form a vector. the mean of each class in the feature space is defined just like how we had discussed in the LDA section.

$$m_i^{\phi} = \frac{1}{n_i}\sum_{l=1}^{n_i}\phi(\mathbf{x}_l^i)$$

Using the equation for $\vec w$ and $m_i^{\phi}$ we can say,

$$\vec w^{\text{T}} m_i^{\phi} = \frac{1}{n_i}\sum_{j=1}^{n}\sum_{l=1}^{n_i}\alpha_jk(\mathbf{x}_j,\mathbf{x}_l^i) = \mathbf{\alpha}^{\text{T}}\mathbf{M}_i$$
where $M_i$ is a vector whose $j^{th}$ element is defined as,

$$(\mathbf{M}_{i})_j = \frac{1}{n_T}\sum_{l=1}^{n_i}k(\mathbf{x}_j,\mathbf{x}_l^i)$$

Here $\mathbf{x}_j$ is the $j^{th}$ data vector among ALL the vectors while $\mathbf{x}_l^i$ is the $l^{th}$ vector in the $i^{th}$ class.

If we define a matrix M as defined below,
$$\begin{align} M & = \sum_{i=1}^k n_i(\mathbf{M}_i-\mathbf{M}_{*})(\mathbf{M}_i-\mathbf{M}_{*})^{\text{T}} \\ \end{align}$$

Where the $j^{th}$ element in the vector $\mathbf{M}_{*}$ is defined as,

$$(\mathbf{M}_{*})_j = \frac{1}{n_T}\sum_{l=1}^{n_T}k(\mathbf{x}_j,\mathbf{x}_l)$$

Then the numerator of the utility function can be written as,
$$\vec w^T S_B^{\phi} \vec w = \alpha^T M \alpha$$

Here $S_B^{\phi}$ is just the between-class covariance matrix defined for the transformed data vectors (in the feature space).

Moving on to the denominator, if we define a matrix $N$ as

$$\begin{align} N & = \sum_{i=1}^k\mathbf{K}_i(\mathbf{I}_{n_i \times n_i}-\mathbf{1}_{n_i})\mathbf{K}_i^{\text{T}} \end{align}$$

with the $(n^{\text{th}},m^{\text{th}})$ component of $\mathbf{K}_i$ defined as $k(\mathbf{x}_n,\mathbf{x}_m^i)$, Also here $1_{n_i}$ is the matrix with all entries $\frac{1}{n_i}$.

Here $\mathbf{K}_i$ is known as the kernel matrix for class $i$.

Then the denominator of the utility function can be written as,
$$\vec w^T S_W^{\phi} \vec w = \alpha^T N \alpha$$

Thus using both of the highlighted equations the utility function in the feature space is,

$$J(\alpha) = \frac{ |\vec w^TS_B^{\phi} \vec w|}{ |\vec w^TS_W^{\phi}\vec w|} = \frac{|\alpha^T M \alpha|}{|\alpha^T N \alpha|}$$

This problem can be solved (just like how LDA was solved in the input space) and the optimal value for $\alpha$ is the leading eigenvector of $N^{-1}M$.

MATLAB Implementation

The above equations can be implemented in MATLAB as follows. Note there is 2 features and we wish to reduce it a single dimension. There are 4 classes ($k = 4$) each class having 13 vectors ($n_i \forall i = 13$) and so the total number of vectors ($n_T = 52$).

• As usual I’ll go ahead and use the Gaussian kernel function. then initialize all the variables.

load TRAINTEST2D

cluster1 = TRAIN{1,6}{1,1}; % Green
cluster2 = TRAIN{1,6}{1,2}; % Blue
cluster3 = TRAIN{1,6}{1,3}; % Red
cluster4 = TRAIN{1,6}{1,4}; % Cyan

nT = 52; % Total number of data Vectors
n1 = 13; % # vectors in cluster 1
n2 = 13; % # vectors in cluster 2
n3 = 13; % # vectors in cluster 3
n4 = 13; % # vectors in cluster 4

% Plot the data before dimensionality reduction
figure(1);
scatter(cluster1(1,:), cluster1(2,:), 'g'); 
hold on;
scatter(cluster2(1,:), cluster2(2,:), 'b');
hold on;
scatter(cluster3(1,:), cluster3(2,:), 'r');
hold on;
scatter(cluster4(1,:), cluster4(2,:), 'c');
hold on;
legend('cluster 1','cluster 2','cluster 3','cluster 4');

• Now I’ll construct the Gram Matrix for all the data vectors,

% CONSTRUCT KERNEL MATRIX
gamma = 10;

% Arrange all the clusters together
X = [cluster1';cluster2'; cluster3'; cluster4'];
K = ones(nT,nT);

% Construct the Gram Matrix 
% K(i,:) - gives the ith vector out of ALL vectors
for i = 1:1:nT
    for j = 1:1:nT
        K(i,j) = kernelGauss(X(i,:),X(j,:),gamma);
    end
end

• Following the equations given above we can construct M,

% Okay K(i,j) = K(x_i,x_j) Now as X(1:13,:), X(14:26,:), X(27:39,:), 
% X(40:52,:) are the 4 clusters respectively. We can calculate M_i as
% the numbers are like that due to how X was constructed.

M_1 = ones(nT,1);
M_2 = ones(nT,1);
M_3 = ones(nT,1);
M_4 = ones(nT,1);
M_star = ones(nT,1);

% sum(K(1:13,j)) is the sum of k(x,x_j) over all x belonging to cluster 1
% sum(K(14:26,j)) is the sum of k(x,x_j) over all x belonging to cluster 2
% sum(K(27:39,j)) is the sum of k(x,x_j) over all x belonging to cluster 3
% sum(K(40:52,j)) is the sum of k(x,x_j) over all x belonging to cluster 4
% sum(K(:,j)) is the sum of k(x,x_j) over ALL x (all clusters)

for j = 1:1:nT
    M_1(j) = (1/nT)*sum(K(1:13,j));
end
for j = 1:1:nT
    M_2(j) = (1/nT)*sum(K(14:26,j));
end
for j = 1:1:nT
    M_3(j) = (1/nT)*sum(K(27:39,j));
end
for j = 1:1:nT
    M_4(j) = (1/nT)*sum(K(40:52,j));
end

for j = 1:1:nT
    M_star(j) = (1/nT)*sum(K(:,j));
end

% Thus we can construct M
M = n1*(M_1-M_star)*(M_1-M_star)' + n2*(M_2-M_star)*(M_2-M_star)' + n3*(M_3-M_star)*(M_3-M_star)' + n4*(M_4-M_star)*(M_4-M_star)';

• For constructing N, we will first separate the gram matrix we had made into the kernel matrices for each of the classes.

% Now we shall construct the kernel matrices for each of the clusters
K_1 = K(:,1:13);
K_2 = K(:,14:26);
K_3 = K(:,27:39);
K_4 = K(:,40:52);

% Thus we can construct N
N = K_1*(eye(n1,n1)-(1/n1)*ones(n1,n1))*K_1' + K_2*(eye(n2,n2)-(1/n2)*ones(n2,n2))*K_2' + K_3*(eye(n3,n3)-(1/n3)*ones(n3,n3))*K_3' + K_4*(eye(n4,n4)-(1/n4)*ones(n4,n4))*K_4';

• Now compute $N^{-1}M$ and its eigenvalues - eigenvectors. Note that as N is usually singular we add a multiple of $I$ for inv(N) to become computable.

% Note that in practice, \mathbf{N} is usually singular and 
% so a multiple of the identity is added to it
N = N + 2*eye(nT,nT);
cond(N) % is N ill conditioned? 

% Now that we have both N and M, we can find all the eigenvectors
% and choose most prominent eigenvectors
P = inv(N)*M;
[V,D] = eig(P);

• Now we can project the data onto a line using the leading eigenvector which has a eigenvalue of 3.07.

% the first eigenvector has the largest corresponding eigenvalue - 3.07
% after which 0.4, 0.28, 0.26 Then it drops down...
alpha = V(:,1);

% the projected points
y = ones(nT,1);

for i = 1:1:nT
    y(i) = alpha'*K(:,i);
end

projCluster1 = y(1:13);
projCluster2 = y(14:26);
projCluster3 = y(27:39);
projCluster4 = y(40:52);

% Plot the projected data
figure(3);
scatter(projCluster1, zeros(13,1), 'g'); 
hold on;
scatter(projCluster2, zeros(13,1), 'b');
hold on;
scatter(projCluster3, zeros(13,1), 'r');
hold on;
scatter(projCluster4, zeros(13,1), 'c');
hold on;
legend('cluster 1','cluster 2','cluster 3','cluster 4');

The plot we get for gamma (of the kernel function) being 10 is,

If we are able to try different values for the kernel function we can surely get great separation of classes.

So in conclusion both methods can be viewed as techniques for linear
dimensionality reduction. However, PCA is unsupervised and depends only on the data vectors (maximize its variance using less features) whereas Fisher linear discriminant (KLDA) also uses class-label information to bring the data vectors closer to other data vectors in the same class while maximizing the variance between classes.

Written with StackEdit.

Equivalent Kernel

The prediction of the output for a input $x_o$ is,
$$y(x_o,m_N) = \vec w_{optimal}^T\vec \phi(x_o)= m_N^T\vec \phi(x_o) = \vec \phi(x_o)^T m_N$$
$$= \vec \phi(x_o)^T\beta S_N\Phi^T\vec t = \beta \sum_{i=1}^n\vec \phi(x_o)^T\Phi^TS_N\vec \phi(x_i)^Tt_i = \sum_{i=1}^nk(x_o, x_i)t_i$$

here the function $k$ is called smoother matrix or the equivalent kernel.

Note that we are substituting the solution for $\vec w$ which is $m_N$ from when we maximize the posterior p.d.f

This effective kernel defines the weights by which the training
set target values are linearly combined in order to make a prediction for a new value of x.

Here scalar $k(x_o, x_i)$ can be thought of as a measure of how important the $i$th input-output vector pair ( the $x_i \rightarrow t_i$ mapping in the Training DataSet) is; for the prediction of the output corresponding to input $x_o$.

For example in the gaussian kernel, the closer $x_o$ is to $x_i$, the larger the magnitude of the kernel function. Thus those corresponding $t_i$ will be weighted higher.

This is visually shown below. when we work with an actual dataset.

You can read more about this in the wiki page for Kernel Smoother.

Kernel Matrix

A Kernel Matrix - $K_{M \times N}$ operates on the training set and produces $M$ predictions, $M$ is the number of input vectors ($x_i^1, x_i^2, x_i^3, ... x_i^M)$, and there are $N$ data samples {$(x_o^1,t_o^1), (x_o^2,t_o^2), (x_o^3,t_o^3), ... (x_o^N,t_o^N)$} in the training dataset.

$$K = \begin{array}{cc} k(x_o^1,x_i^1) & k(x_o^1,x_i^2)& ... & k(x_o^1,x_i^N) \\ k(x_o^2,x_i^1)& k(x_o^2,x_i^2) & ... & k(x_o^2,x_i^N)\\ ... & ... & ... & ... \\ k(x_o^M,x_i^1)& ... & ... & k(x_o^M,x_i^N)\\ \end{array}$$

$$\vec t_{N \times 1} = [t_o^1, t_o^2, t_o^3, ... t_o^N]^T$$

Here the predictions for the input vectors will be a vector - $\vec p$,

$$\vec p = K\vec t$$

If we wish to form a “model” of sorts (which gives us predictions) instead of laboriously computing the kernels every time we wish to make a prediction. We can form the kernel and thus for any new input we can find the closest $x_o$. In order to save memory space and even ensure accuracy we can set a resolution $r$ and operate on normalized data.

So if we consider a resolution of 0.01 and the normalization range to be [0,1]. Then the input vectors ($x_i^1, x_i^2, x_i^3, ... x_i^M)$ will be $(0,0) , (0,0.01), ... (0.01,0) , (0.01,0.01), .... .... (1,1)$ if we consider the input to be 2 dimensional.

So even though it will take a long time to compute the entire matrix and even experiment with various kernel parameter values to finally obtain the kernel matrix. But once we have it we can quickly finish predicting for new inputs.

Now I’ll use the above concepts to form such a model using this dataset. Find the whole code here, please ensure you have kernel Functions defined in separate files and added to MATLAB’s path.

% The exponential kernel is closely related to the Gaussian kernel, with only the 
% square of the norm left out. It is also a radial basis function kernel.
function [res] = kernelExp(X1, X2, c)
    res = exp( -(norm(X1-X2)) * c );

• First load the data into the workspace,

load('ASSIGNMENT1.mat');

• Separate DATA into the training dataset and testing dataset.

% Training DataSet - 7,500 data samples
inputTrain = DATA(1:7500,1:2); % each row is a 1*2 input vector
outputTrain = DATA(1:7500,3:4); % each row is a 1*2 target vector

% Testing DataSet
inputTest = DATA(7501:10000,1:2);
outputTest = DATA(7501:10000,3:4);

• Now depending on how the input is normalized (in our case it’s 0 to 1), get the parameters initialized. The larger the resolution the longer the forming of the model will take, but the model will serve to predict a larger number of inputs with great accuracy (we won’t have to round off to the closest). I am preallocating the modelInputsmatrix (as when an array is growing by assignment or concatenation it affects the code performance).

r = 0.1; % Resolution
% The number of X_o's (as input is 1*2) and normalized to [0,1]
M = int8( ((1/r)+1)*((1/r)+1) ); 
N = 7500; % 75% of the total 10,000 samples are used as the training set.
modelInputs = zeros(M,2);
NormalisedInputs = 0:r:1;

• Construct the set of input vectors using NormalisedInputs, by taking all pairwise combinations of it (as we need a 1*2 vector).

for i = 1:1:length(NormalisedInputs)
    for j = 1:1:length(NormalisedInputs)
        %modelInputs = [modelInputs; [NormalisedInputs(i),NormalisedInputs(j)]];
        modelInputs(j+length(NormalisedInputs)*(i-1),1:2) = [NormalisedInputs(i),NormalisedInputs(j)];
    end
end

• For a particular value of gamma now we can construct the matrix,

K = zeros(M,N);
gamma = 0.1;

for i = 1:1:M
    for j = 1:1:N
        K(i,j) = kernelGauss(modelInputs(i,1:2),X_i(j,1:2),gamma);
    end
end

• Using the Gaussian basis function and the Sigmoidal basis function, we can compute the kernel matrix (smoothing matrix).

They are the following,
$$\Large k_{\gamma}(\vec x_o,\vec x_i) = e^{- \gamma (|| \vec x_o - \vec x_i||^2) }$$
$$\Large k_{\alpha, c}(\vec x_o,\vec x_i) = \tanh (\alpha \vec x_o^T \vec x_i + c)$$

K_gauss = zeros(M,N);
K_sigmoid = zeros(M,N);
gamma = 10;
alpha = 1;
c = 1;

for i = 1:1:M
    for j = 1:1:N
        K_gauss(i,j) = kernelGauss(modelInputs(i,1:2),X_i(j,1:2),gamma);
        K_sigmoid(i,j) = kernelHyperTangent(modelInputs(i,1:2),X_i(j,1:2),alpha,c);
    end
end

• At this point I would like to show visually how the weighing of the training set works, as you can see each row of the $K$ matrix corresponds to a single input while the columns are for each of the training samples. So if you run the following code,

Note that I tried for about 2 days to find the MATLAB syntax to represent the X axis as the first input parameter, Y axis as the second input parameter. So when we take say resolution of 0.1. then that would be the scale of the X and Y axis. Now using Z axis we can represent the weight we attribute to that input. So it would be a scatter plot with (x,y,z) where (x,y) is the input and z is its corresponding weight. In this way it would be clearly seen as a hill how the inputs similar to the training input under consideration gets high weightage.

The following code displays the weights corresponding to all 121 input vectors from (0,0) to (1,1) (when we took resolution = 0.1) for the Pth training sample. Here M is 121 cause that many $1\times2$ input vectors exist from 0 to 1 with a resolution of 0.1

plot(K_gauss(1:M,P));

When P = 98 (the index of the training sample is 98) with training input - (0.396, 0.495), we get the following plot.

See how we have 5 significant peaks in the weights? if we had plotted it in 3d we could have seen all these peaks were actually close together.
The largest peak is at index 50 which happens to be the input vector (0.4,0.5). Note how similar it is to the training input we are comparing with.
The next 2 peaks are at 39 and 61, which are (0.3,0.5) and (0.5,0.5). So we see a unit resolution change also counts as close.
The final 2 are 28 and 72 which are (0.2,0.5) and (0.6,0.5).

So this shows us how $k(x_o, x_i)$ is larger the closer $x_o$ is to $x_i$.

Now the question arises how much larger should the weights be, if they are close, rather how should the weights vary as the distance between them - $(x_o, x_i)$ changes. This is controlled using the $\gamma$ parameter, I used a value of 50 and obtained the above graph. Note that a larger gamma, means that more emphasis on the closeness. So for $\gamma = 500$, you will get a single peak at the the 50th index.

This much must have got you a great idea on how this gaussian kernel works as an instance based learner.

Ok so for a particular training sample (i.e. $98^{th}$) we plotted the weight values (corresponding to each of the input vectors - 121 of them) and showed how the input vectors which are “closer” to that particular training sample’s input is given a larger weight.

Cool so to predict the output for a particular input vector - we use all the weights and weigh ALL the training sample’s outputs.

But first we need to normalize the weights to 1 as we use these weights to only contribute the relative importance of the various training set’s outputs.

% This vector's ith element is the sum of weights used to predict the 
% output for the ith input vector.
sumOfWeights = zeros([M,1]);
for i = 1:1:M
    sumOfWeights(i) = sum(K_gauss(i,1:N));
end

predictions = K_gauss * outputTrain;
% Normalize it
predictions(1:M,1) = predictions(1:M,1)./ sumOfWeights;
predictions(1:M,2) = predictions(1:M,2)./ sumOfWeights;

Moving on to testing this model, we can calculate the MSE when we compare the results of the prediction to the actual outputs.

% Now we are going to test this model - calculate MSE error
L = 10000 - N; % remmaining samples in the data set are used to test
mse = 0;

for i = 1:1:L
    % in this iteration we wish to find the output for inputTest(i,1:2)
    modelInput = inputTest(i,1:2);

    % First we need to find which index this input is closest to
    % in modelInputs
    X1 = round(modelInput(1), numberOfDecimals);
    X2 = round(modelInput(2), numberOfDecimals);

    index1 = find(modelInputs(1:M,1) == X1);
    index2 = find(modelInputs(1:M,2) == X2);
    index = intersect(index1,index2);

    % so the corresponding prediction will be
    modelOutput = predictions(index,1:2);
    actualOuput = outputTest(i,1:2);

    squaredError = (modelOutput - actualOuput)*(modelOutput - actualOuput)';
    mse = mse + squaredError;
end

% as its the Mean we gotta divide with number of samples.
mse = mse/L;

With $\gamma = 50$ using the gaussian kernel we get a mean square error of $1.1987$ which is reasonable.

Written with StackEdit.

Bayesian Technique

Here we will explore the relationship between maximizing likelihood p.d.f - $p(\vec t|\vec w)$, maximizing posterior p.d.f - $p(\vec w|\vec t)$, minimization of the sum-of-squares error function - $E_D(\vec w)$ and the regularization technique.

When we maximise the posterior probability density function w.r.t the parameter vector using the “bayesian technique” - we need both the likelihood function and the prior. (The denominator in the bayes theorem is just a normalization constant so that doesn’t really matter.)

$$p(\vec w|\vec t) \propto p(\vec t|\vec w)p(\vec w)$$

The model

We have a set of inputs, $\vec x= [x_1, . . . , x_n]^T$ with corresponding target values $\vec t =[t_1, .... t_n]^T$.

We assume that there exists some deterministic function $y$ such that we can model the relationship between these two as the sum of $y(x,w)$ with additive gaussian noise,

$$t_i = y(x_i,\vec w) + \mathcal{N}(0,\beta^{-1})$$

$β$ is the precision (inverse variance) of the additive univariate Gaussian noise.

We define $y$ as the linear combination of basis functions,
$$y(x_i,\vec w) = w_1\phi_1(x_i)+w_2\phi_2(x_i)+...+w_p\phi_p(x_i) = \vec w^T\vec \phi(x_i)$$
We define the parameter vector as $\vec w_{p \times 1}=[w_1,w_2,...,w_p]^T$ and basis vector as $\vec \phi(x_i)=[\phi_1(x_i),\phi_2(x_i),...,\phi_p(x_i)]^T$.

This parameter vector is very important as the posterior p.d.f is the updated probability of $\vec w$ given some training data. which is found from the prior data of $\vec w$. While the likelihood p.d.f of getting that training data given $\vec w$.

We usually choose $\phi_1(x)=1$ because we need a bias term in the model. (to control the extent of the shift in $y$ itself - check this answer out)

For the data set as a whole we can write the set of model outputs as a vector $\vec y_{n \times 1}$,

$$\vec y(\vec x,\vec w) = \Phi\vec w$$

Here the basis matrix $\Phi_{n \times p}$ is a function of $\vec x$ and is defined with its $i$th-row being = $[\phi_1(x_i),\phi_2(x_i),...,\phi_p(x_i)]$ for $n$ such rows.

Likelihood function

We assume that these data points $(x_i,t_i)$ are drawn independently from the distribution we would have to multiply the individual data point’s p.d.f - which are gaussian.

$$p(\vec t|\vec x,\vec w, β) =\prod^n_{i=1}\mathcal{N}(t_i|\vec w^T\vec \phi(x_i), β^{-1})$$

Note that the $i$th data points p.d.f is centered around $\vec w^T\vec \phi(x_i)$ as the mean.

Does the product of $n$ univariate gaussians forms a multivariate distribution in {$t_i$}?? I say this because we choose a gaussian prior, thus the likelihood should also be gaussian right?

Prior

We choose the corresponding conjugate prior, as we have a likelihood function which is the exponential of a quadratic function of $\vec w$.

No clue why but for now for this to make sense let’s say that the likelihood function is also gaussian - product of all those gaussians.

Thus the prior p.d.f is a normal distribution - $\mathcal{N}(m_o,S_o)$

Posterior

The posterior p.d.f is a $\mathcal{N}(m_N,S_N)$ (as we choose a conjugate prior)

After solving for $m_N$ and $S_N$ we get,

(The complete derivation is available in Bishop - (2.116)) - coming soon

$$m_N = S_N(S^{-1}_0 m_0 + βΦ^T\vec t)$$
$$S^{-1}_N = S^{-1}_0 + βΦ^TΦ$$

The sizes are,
The mean vectors, $m_N$ and $m_o$ are both $p \times 1$ and they can be thought of as the optimal parameter vector and pseudo observations respectively.
The covariance matrices, $S_N$ and $S_o$ are both $p \times p$

We shall consider a particular form of Gaussian prior in order to simplify the treatment. Specifically, we assume a zero-mean isotropic Gaussian governed by a single precision parameter $α$,
$$\mathcal{N}(\vec 0,α^{-1}I)$$
So we basically take $m_o = \vec 0$ and $S_o = α^{-1}I_{p \times p}$

Thus if we use this prior we can simplify the mean vector and the covariance matrix of posterior p.d.f to,

$$m_N = βS_NΦ^T\vec t$$
$$S^{-1}_N = αI + βΦ^TΦ$$

Now if we take log of the posterior pdf - $\mathcal{N}(m_N,S_N)$ , in order to maximize it with respect to w, we find that what we obtain is equivalent to the minimization of the sum-of-squares error function with the addition of a quadratic regularization term, corresponding to $λ = \frac{α}{β}$.

$$\ln p(\vec w|\vec t) = \frac{-\beta}{2}\sum_{i=1}^N(t_i - \vec w^T\vec \phi(x_i))^2 - \frac{α}{β}\vec w^T\vec w = -βE_D(\vec w)-\frac{α}{β}\vec w^T\vec w$$

Thus we conclude that while maximising likelihood function is equivalent to the minimization of the sum-of-squares error function, maximising posterior p.d.f is equivalent to the regularization technique.

The regularization technique is used to control the over-fitting phenomenon by adding a penalty term to the error function in order to discourage the coefficients from reaching large values.

This penalty term arises naturally when we maximize posterior p.d.f w.r.t $\vec w$

Here the minimization of the sum-of-squares error function - $E_D(\vec w)$ is also same as Maximization of the likelihood p.d.f. Taking log of $\ln p(t|\vec w, β)$ we get,

$$\ln p(t|w, β) = \frac{N}{2}\ln \beta - \frac{N}{2}\ln(2π) - \beta E_D(\vec w)$$

thus maximizing likelihood is equal to maximizing $- E_D(\vec w)$ (rest are all constants w.r.t $\vec w$)

Parameter Estimation

Objective

Use System Identification Toolbox of MATLAB to estimate the transfer function of a system - given sampled data of the system’s input and output.
Then use the minimum mean square error (MMSE) estimator method to alternatively estimate the model of the system.

Theory

Mathematical modelings importance in various fields cannot be understated. It is used to describe a system using mathematical concepts and language. When we have an interest in performing tests on a complicated system which is difficult to physically obtain, expensive, and sensitive to failure. Then it is safer and cheaper to perform the same tests on the model using computer simulations rather than carry out repetitive experimentations and observations on the real system.

In making these “Models” we have various types,

Theoretical Models

These models are obtained from fundamental principles, such as the laws of conservation of mass, energy, and momentum along with other chemical principles such as chemical reaction kinetics and thermodynamic equilibrium, etc.

Empirical models

These models are based on experimental plant data. These models are developed using data fitting techniques such as linear and nonlinear regression.

We will be dealing with empirical models in the remainder of this experiment, models obtained exclusively from experimental plant data are also known as black-box models . Such models do not provide detailed description of the underlying physics of the process. However, they do provide a description of the dynamic relationship between inputs and outputs. Thus they are sometimes more adequate for control design and implementation.

Thus we come to the issue of Parameter Estimation, given a set of input output data pairs we need to estimate the values for the coefficients used in the model such that the model best fits the experimental evidence.

Let the relationship between input and output of the system be modeled as the difference equation given below,

$$y_{k+1} = a.y_k + b.u_k$$

where $y_k$ is the $k^{th}$ output sample and $u_k$ is the $k^{th}$ input sample.

So for $n$ samples of input-output data pairs, we can write the set of simultaneous equations in matrix form in order to form a concise notation.

Consider a $(n\times1)$ vector $Y$ having elements from $y_1$ to $y_n$,
and a $(n\times2)$ matrix $\Phi$ having elements from $y_0$ to $y_{n-1}$ as the first column and $u_0$ to $u_{n-1}$ as the second column,
and finally a $(2\times1)$ vector $\theta$ which holds the two parameters $a$ and $b$.

MMSE Derivation

So the $n$ equations formed from entering the n data pairs into the difference equation given above can be written in matrix form as follows,

$$\Large Y = \Phi \theta$$

even though we write a equals sign, it is usually not possible to find a $\theta$ which perfectly models the system. we just aim for a $\theta$ such that the model most faithfully fits the given data ($n$ samples). Note that as $lim_{n \to \infty} Y = \Phi \theta$, (infinite number of data points) until then R.H.S and L.H.S continues to approach each other.

From the minimum mean square error (MMSE) estimator method, we get the following formula to calculate this $\theta$,

$$\Large \theta = (\Phi^T.\Phi)^{-1}.\Phi^T.Y$$

I will now detail the proof for the above equation which is used to estimate $\theta$ (i.e. the values of $a$ and $b$).

As I said before we need to find a $\theta$ such that the model most faithfully fits the given data, how can we make this statement into a relationship involving the parameters of interest?

So “most faithfully” can be thought of to mean, we want to minimize the error. Let us then define the error of the model as follows,

$$E = O_{SystemsOutput} - O_{ModelsOutput} = Y - \Phi \theta$$

Now we don’t want to minimize this error - $E$ for all data samples. Because we care about the deviation of $O_{Models-Output}$ from $O_{Systems-Output}$, we do not care about the direction of this deviation. This is because positive and negative error of subsequent samples can cancel each other out. Thus we should take the $|E|$ - so that we have the absolute value of Error. We would rather define the error as the value $P$ (also known as Mean Square Error) rather than just $|E|$.

$$P = \frac{E^T E}{n} =\frac{1}{n} (Y_{n\times1} - \Phi_{n\times2} \theta_{2\times1})^T(Y_{n\times1} - \Phi_{n\times2} \theta_{2\times1})$$

The reason why we prefer MSE (squaring) to taking modulus when we want to quantify this error, is due to the fact that while both eliminate the primary problem of “positive and negative error”; squaring makes the algebra much easier to work with and offers properties that the absolute method does not. Additionally the squared difference has specific mathematical properties; such that it’s continuously differentiable (which nice when you want to minimize it - like now). There are additional reasons which can be found here.

Okay getting back to our hunt for the best fit for the given data, we need to minimize $P$ with respect to $\theta$, for that we differentiate $P$ with respect to $\theta$ and equate that to zero.

Lets expand $P$,
$$P = \frac{1}{n} (Y - \Phi \theta)^T(Y - \Phi \theta)$$
The transpose respects addition,
$$P = \frac{1}{n} (Y^T - (\Phi \theta)^T)(Y - \Phi \theta)$$
Opening the brackets, note that $(\Phi \theta)^T = \theta^T\Phi^T$
$$P = \frac{1}{n} (Y^TY - Y^T\Phi\theta - \theta^T\Phi^{T}Y + \theta^T\Phi^T \Phi \theta)_{1\times1}$$

Vector Calculus Side Note
$$\frac{d(x^Ta)}{dx} = \frac{d(a^Tx)}{dx} = a^T$$
$$\frac{d(x^TAx)}{dx} = x^T(A^T + A)$$
Here both $x$ and $a$ are vectors while $A$ is a matrix.

---

You can skip the following proof

Now the formulas given above can be understood by the fact that vectors are just used to hold multiple elements together - for the sole reason of applying operators onto them in bulk.

So in the first formula both $a$ and $x$ are $n\times 1$ vectors. What both $x^Ta$ and $xa^T$ will be is a scalar value which is the dot product of both vectors.
$$x = [x_1, x_2, ... x_n]^T$$
$$a = [a_1, a_2, ... a_n]^T$$
$$k = x^Ta = xa^T = [x_1a_1 +x_2a_2 + ... + x_na_n]_{1\times 1}$$
Differentiating a scalar w.r.t a vector is defined as a new vector with its $i$th element being the differentiation of that scalar with the $i$th element of the vector.

Here I have to add that this topic divides the academia.
Two competing notational conventions split the field of matrix calculus into two separate groups. The two groups can be distinguished by whether they write the derivative of a scalar with respect to a vector as a column vector or a row vector.
$$\frac{d(x^Ta)}{dx} = \frac{d(a^Tx)}{dx} = [\frac{dk}{dx_1}, \frac{dk}{dx_2}, ... \frac{dk}{dx_n}]$$
$$=[a_1, a_2, ... a_n ] = a^T$$

For the second formula, the proof is straightforward given the above result. Here $A$ is a $n \times n$ matrix. Apply product rule, note that the terms in the bracket are considered as constants when differentiating.
$$\frac{d(x^TAx)}{dx} = \frac{d((x^TA)x)}{dx} + \frac{d(x^T(Ax))}{dx}$$

From the first formula - when comparing, $(x^TA)_{1 \times n}$ can be thought of as $a^T$ while $(Ax)_{n \times 1}$ as $a$. So,
$$= x^TA + (Ax)^T = x^TA + x^TA^T = x^T(A^T + A)$$
$$Q.E.D$$

---

When we partially differentiate P w.r.t $\theta$, only $\theta$ is a variable while everything else is considered to be a constant.

Partially differentiating w.r.t $\theta$ and equating to 0 we get the value of $\theta$ for which the Mean Square Error is minimum,
$$\frac{\partial P}{\partial \theta} = \frac{1}{n} ( [0]_{1\times2} - Y^T \Phi - (\Phi^T Y)^T + \theta^T((\Phi^T \Phi)^T + \Phi^T \Phi) ) = [0]_{1\times2}$$

Here $\Phi^T \Phi$ is a symmetric matrix thus,
$$\frac{\partial P}{\partial \theta} = \frac{1}{n} ( [0]_{1\times2} - Y^T \Phi - (\Phi^T Y)^T + 2\theta^T( \Phi^T \Phi) ) = [0]_{1\times2}$$

Moving all terms without $\theta$ to one side,
$$2\theta^T( \Phi^T \Phi) = 2 Y^T\Phi$$

So post-multiplying both sides with $( \Phi^T \Phi)^{-1}$,
$$\theta^T = Y^T\Phi( \Phi^T \Phi)^{-1}$$

Taking transpose on both sides, note that $(( \Phi^T \Phi)^{-1})^T = (( \Phi^T \Phi)^{T})^{-1}$ because ${(A^{-1})}^T= {(A^T)}^{-1}$ and $(AB)^T = B^TA^T$.

Thus we get our $\theta$,
$$\theta = (\Phi^T.\Phi)^{-1}.\Phi^T.Y$$

Simulation

I will first generate “data” using the simulink library in MATLAB. There we can construct a model of a system and run simulations on it.

Here as you can see I am applying a step input to a first order system, and writing both the input (Step) and output to the file try.mat, now looking at the settings of the step function and the ‘write to file’ block you can see I have chosen a Sampling Rate of 0.1 secs and initial time is 0 secs.

So we store time, input, and output into a variable - data and this file is called try.mat. Note that as we have 0.1 sec sample time we sample every 0.1 secs for 30 secs. Therefore we will have 301 samples ($\frac{30}{0.1}$ and an initial 0 data sample).

Here the transfer function of the system used to generate the data is ,
$$G_p = \frac{5}{3s+1}$$

You might get a error - Error encountered in block ‘exp1/To File’ - Error opening or closing file ‘try.mat’ . The error will open the Diagnostic Viewer as shown below,

We can fix this error by opening MATLAB as administrator, this error is caused due to MATLAB not having the permission to write to disk. Open MATLAB as shown in the image below,

On opening the Scope block we will see the input and output superimposed on a graph with amplitude as a function of time. The yellow line is the unit step input and the purple line is the systems output.

In order to do that we need to load this generated data into the MATLAB’s workspace. We use the load try.mat command to achieve this, once we do that the data variable holding the time, input, and output as rows will be added to the MATLAB workspace. As you can see below the first row is the time sampled every 0.1 secs, the second row is the step input which starts at 1 sec and finally the system’s output is the third row in the data variable. Note that I printed data' for better readability that is why the rows I refer to are shown as columns in the image below.

This data variable can now be separated into its 3 components using the commands as shown below,

We will be now using the t, u, and y variables we have made in order to reverse engineer the above system transfer function.

The Minimum Mean Square Error (MMSE) estimator method

The first method is using the equation we just derived which models the system using a first order difference equation, and then minimizes the Mean Square Error between the actual system’s output and the model’s predicted output.

Now that we have the time, input, and output variables in the workspace, we can construct $Y$ and $\Phi$ using the commands as shown below,

Note that in MATLAB element’s indexes range from 1 to n rather than 0 to n like more programming languages.

As I had explained in the proof, I am taking the all but the first data element when I construct $Y$, and all but the last element when i construct $\Phi$.

Then I apply the formula to find the $\theta$ for which the MSE is minimum. This comes out to be for values of $a$ = 0.9618 and $b$ = 0.1950.

Let us now verify how well these parameters fit the actual system, for this I am using MATLAB’s c2d function to take my continuous time i/o system and find its equivalent discrete time i/o system with a sampling time of 0.1 secs (Just like the data we generated).

Note that the $dG$ shown in the image is a approximation of $cG$, and with smaller sampling time this approximation eventually tends to become the continuous-time model itself!

For better understanding I will elaborate on how the $dG$ tells us what $a$ and $b$ are. MATLAB tells us,

$$dG = \frac{Y(z)}{U(z)} = \frac{0.1639}{z - 0.9672}$$

Cross multiplying, we get
$$z Y(z) - 0.9672Y(z) = 0.1639 U(z)$$

By the Time shifting property of Z transforms i.e. $z^{-k}X(z) = x[n-k]$,
$$Y(z+1) = 0.9672Y(z) + 0.1639U(z)$$

Doesn’t the above equation look familiar? yes it is of the same structure as the difference equation we used to model this system. Just as first order systems are modeled using $\frac{k}{\tau s + 1}$, we can use the difference equation of this form to model the equivalent discrete time model system.

Thus the best fitting which MATLAB performs using its “inbuilt algorithms” (which are based on the same mathematics we used, but has multiple small modifications to improve its speed and accuracy) is $a$ = 0.9672 and $b$ = 0.1639.

So an error of $\delta a$ = 0.0054 and $\delta b$ = -0.0311 seems acceptable.

System Identification Toolbox

We can also use the System Identification Toolbox that comes with MATLAB to estimate the same parameters. This method has the additional benefits in the fact that they have added options for all algorithms to estimate the parameters, runs extremely fast, graphs the model formed etc.

We can use ident to open the System Identification App, we will see a window as shown below open up,

We can now import the data using which MATLAB will estimate the parameters,

Note that we need to keep the sampling time as 0.1 sec and starting time as 0 secs just as we defined it in the simulink model. Now click import to get the data into ident. The time domain data will look like the below image,

---

Now using Estimate -> Transfer Function Models, we are prompted for number of poles and zeros which we need to enter.

After which they estimate using Instrument Variable approach (which is the default) as shown below,

Now as you can see the Transfer function has been estimated with a 100% fit.

This is absolutely correct as $\frac{5}{3}$ is 1.666 and $\frac{1}{3}$ is 0.333.

We can even use the “Show in LTI” option as shown in image above to see the out of a LTI system modeled on a $G(s) = \frac{\frac{5}{3}}{s+ \frac{1}{3}}$.

---

There is also a different MATLAB option to do the estimation, that is using
Estimate -> Process Models,

Just untick the ‘Delay’ checkbox, uncheck ‘Zero’ checkbox - as there are no zero’s, choose the correct number of poles - 1 from the dropdown menu. Then click on estimate - We get the results quite fast,

As we can see the transfer function has been recovered from the time domain data here, using the $\frac{k}{Tp_1 s + 1}$ model. Where $k$ is 5 and $Tp_1$ is 3.

Conclusion

In conclusion we have used the experimental plant data (Which we generated using simulink) to estimate black-box models . We understood the mathematics behind this generation and also employed more suitable methods available in MATLAB in the process of generating said models.

Addendum

We can show very clearly the benefits of empirical models when we model the FET amplifier circuit - voltage divider bias configuration.

Now before we begin, let’s remember that a DC gate-to-source voltage controls the level of DC drain current through a relationship known as Shockley’s equation.
$$I_D = I_{DSS}(1-\frac{V_{GS}}{V_P})^2$$

So we can express the change in drain current that will result from a change in gate-to-source voltage by using the transconductance factor $g_m$.
$$\Delta I_D = g_m \Delta V_{DSS}$$

Now substituting the ac equivalent model for the JFET, we get the following circuit,

From the above circuit we find the input impedance $Z_i$ and output impedance $Z_o$ as,
$$Z_i = R_1||R_2$$
Setting $V_i$ = 0 Volts sets $V_{GS}$ and $g_m V_{GS}$ to zero, also
$$Z_o = r_d||R_D$$

So the Gain of the amplifier $A_v$ is,
$$A_v = \frac{V_o}{V_i} = \frac{-g_m V_{GS}( r_d||R_D)}{V_{GS}} = -g_m V_{GS}\frac{r_dR_D}{r_d+R_D}$$

The analysis thus far has been limited to a particular frequency. For the amplifier, it was a frequency that normally permitted ignoring the effects of the capacitive elements, reducing the analysis to one that included only resistive elements and sources of the independent and controlled variety.

The junction capacitance (input capacitance) $C_m$ will be virtually increased due to Miller’s effect.
$$C_m = C(1+A_v)$$
where $-A_v$ is the gain of the amplifier and $C$ is the feedback capacitance.
$$-A_v = \frac{V_o}{V_i'} = \frac{V_o}{\frac{V_i}{1+j\omega R_a C_m}}$$
So now the Transfer function is given by,
$$G(s) = \frac{V_o}{V_i} = \frac{-A_v}{1+j\omega R_a C_m} = \frac{-g_m V_{GS}}{1+s R_a C_m}\frac{r_dR_D}{r_d+R_D}$$

Now just like how we obtained the above transfer function by first principle we can use the parameter estimation methods detailed in the sections above to obtain the same transfer function.

Thus knowing that FET amplifier has a first order transfer function characteristics lets us model it using $\frac{k}{\tau s + 1}$.

Note that this is not a complete black box model, it is infact a grey box model as we know the structure (First order).

Written with StackEdit.