Linear Systems Review
Daniel Kersten

Vector operations

Dimension of a vector.

You can get the dimensionality of a vector using  Dimensions[], or Length[].

v = {2.1, 3, -0.45, 4.9};

Dimensions[v]

Dimensions[], will give you the dimensions of a matrix, while  Length[] tells you the number of elements in the list. For example,

M = {{2,4,2}, {1,6,4}};

Length[M]

Try comparing  Length[M] with Dimensions[M].

Transpose of a vector.

The transpose of a column vector is just the same vector arranged in a row. However, because of the way Mathematica uses lists to represent vectors you don't have to distinguish between row and column vectors. In standard math notation, transpose of a vector x, is often written xT. You can see a vector in column form by typing v//MatrixForm, or:

MatrixForm[v]

Vector addition is accomplished by simply adding the components of each vector to make a new vector. Note that the vectors all have the same dimension.

a = {3,1,2};
b = {2,4,8};
c = a + b

Vectors can be multiplied by a constant. We saw an example of this earlier.

2 a

Metric length of a vector.

It is unfortunate terminology, but Length[] does NOT give you the metrical length of the vector. In order to get the length of a vector, you calculate the Euclidean distance from the origin to the end-point of the vector. We get this by squaring each component, adding up the squares, and taking the square root. First, we will do this using the Apply[] function, where the Plus operation is applied to all the elements of the list. Note that the operation of exponentiation is "listable", that is it is applied to each element of the vector:

a^2

What is a a ?

N[Sqrt[Apply[Plus, a^2]]]

If you wish, you can define your own function to apply to the list. What we have just calculated is the square root of the dot product or inner product of a with itself. The length of a vector a is often written as |a| in standard math notation. In the next section, we use the inner or dot product to calculate the metric length of a vector.

Inner product. To calculate the inner product of two vectors, you multiply the corresponding components and add them up:

u = {u1,u2,u3,u4};
v = {v1,v2,v3,v4};
u.v

The inner product is also called the dot product. Later we will see what is meant by outer product. The inner product between two vectors a and b is written either as:
                                                   a.b or [a,b], or aTb
Mathematica  uses the dot notation.

One use of the inner product is to calculate the length of a vector. a.a  is just the sum of the squares of the elements of a, so gives us another way of calculating the length of a vector.

N[Sqrt[a.a]]

Let's define a function that will return the length of a vector, x:

Vectorlength[x_] := N[Sqrt[x.x]]

Projection. The dot product, a.b, is equal to:

                                                   |a| |b| cos(angle between a and b)

In problem set 1, you calculate the output of a linear neuron model as the dot product between an input vector and a weight vector. Both the weight and input lists can be thought of as vectors in an n-dimensional space. Suppose the weight vector has unit length.  Recall that you can normalize any vector to unit length by dividing by its length:

v = v/Sqrt[v.v] ;

Geometrically, we can think of the output of a neuron as the projection of the activity of the neuron input activity vector onto the weight vector direction. Suppose the input vector is already perpendicular to the weight vector, then the output of the neuron is zero, because the cosine of 90 degrees is zero. As you found or will find with the cross-correlator of Problem Set 1, the further the input pattern is away from the weight vector, as measured by the cosine between them, the poorer the "match" between input and weight vectors, and the lower the response.
Here are three lines of code that calculate the two-dimensional vector z in the direction of w,  with a length determined by "how much of x projects in the w direction":

x = {1,2};
w = N[{2/Sqrt[5],1/Sqrt[5]}];
z = (x.w) w ;

Show[   Graphics[{Line[{{0,0}, x}],  
        	Line[{{0,0}, z}],
        	{Dashing[{0.03,0.03}], Line[{x, z}] },
       		Text[FontForm["w",{"Helvetica-Bold",18}], w, {0,-1}],
        	Text[FontForm["x",{"Helvetica-Bold",18}], x, {-2,0}],
        	Text[FontForm["z",{"Helvetica-Bold",18}], z, {0,-1}],
        	{AbsoluteThickness[3], Line[{{0,0}, w}] }
        }], 
        Axes->True, AspectRatio->1
        
];

Angle between two vectors and orthogonality: Similarity measure between patterns

Often we will want some measure of the similarity between two patterns of neural firings. As we have just seen, one measure of comparison is the degree to which the two state vectors point in the same direction. The cosine of the angle between two vectors is one possible measure:

cosine[x_,y_] := x.y/(Vectorlength[x] Vectorlength[y])
cosine[a,b]

Note that if two vectors point in the same direction, the cosine of the angle between them is 1:

a = {2,1,3,6};
b = {6, 3, 9, 18};
cosine[a,b]

Try verifying that w and z from the previous section point in the same direction.
If two vectors point in the opposite directions, the cosine of the angle between them is -1:

a = {-2,-1,-3,-6};
b = {6, 3, 9, 18};
cosine[a,b]

Two vectors may point in the same direction, but could be quite different because they have different lengths. Another measure of similarity is the length of the difference between two vectors:

Vectorlength[a - b]

Orthogonality. The case where vectors are at right angles to each other is an important special case that is worth spending a little time on. Consider an 8-dimensional space. One  very familiar set of orthogonal vectors is the following:

u1 = {1,0,0,0,0,0,0,0};
u2 = {0,1,0,0,0,0,0,0};
u3 = {0,0,1,0,0,0,0,0};
u4 = {0,0,0,1,0,0,0,0};
u5 = {0,0,0,0,1,0,0,0};
u6 = {0,0,0,0,0,1,0,0};
u7 = {0,0,0,0,0,0,1,0};
u8 = {0,0,0,0,0,0,0,1};

Each vector has unit length, and it is easy to see just by inspection that the inner product between any two is zero. On the other hand, here is another set of 8 vectors in 8-space for which it is not immediately obvious that they are all orthogonal. These vectors are called Walsh functions:

v1 = {1, 1, 1, 1, 1, 1, 1, 1};
v2 = {1,-1,-1, 1, 1,-1,-1, 1};
v3 = {1, 1,-1,-1,-1,-1, 1, 1};
v4 = {1,-1, 1,-1,-1, 1,-1, 1};
v5 = {1, 1, 1, 1,-1,-1,-1,-1};
v6 = {1,-1,-1, 1,-1, 1, 1,-1};
v7 = {1, 1,-1,-1, 1, 1,-1,-1};
v8 = {1,-1, 1,-1, 1,-1, 1,-1};

You can calculate the inner products between any two, and you will find out that they are all zero. Note that with the first set of vectors, {ui}, you can tell which vector it is just by looking for where the 1 is. For the second set, {vi}, you can't tell by looking at just one component. For example, the first component of all of the Walsh functions has a 1. You have to look at the pattern to tell which Walsh function you are looking at. 
Suppose for the moment that we want to assign meaning to each of the patterns--each pattern is a code for some thing, like "grandma Tompkins", "grandma Wilke", and so forth. If we use the u's, then we could look for the one neuron that lights up to find out which grandma it is representing--then neuron activity represented, for example, by the third element of the pattern could mean "grandma Wilke". This strategy wouldn't work if we encoded a collection of grandmothers using the v's. The v's give us a simple example of what is sometimes referred to as a distributed code. The w's are examples of a grandmother cell code. The reason for this obscure terminology can be traced to earlier debates on whether there may be single cells in the brain whose firing uniquely determines the recognition of one's grandmother.

Orthonormality. The Walsh set is orthogonal, but they are not of unit length. We have already seen some of the advantages of working with unit length vectors. The general issue of normalization comes up all the time in neural networks both in terms of limiting overall neural activity, and limiting synaptic weights. So it is sometimes convenient to normalize an orthogonal set, producing what is known as an orthonormal set of vectors:

w1 = v1/Vectorlength[v1];
w2 = v2/Vectorlength[v2];
w3 = v3/Vectorlength[v3];
w4 = v4/Vectorlength[v4];
w5 = v5/Vectorlength[v5];
w6 = v6/Vectorlength[v6];
w7 = v7/Vectorlength[v7];
w8 = v8/Vectorlength[v8];

Vector representations, linear algebra

The issue of how information is to be represented is fundamental in the information sciences generally, as well as for neural network theory. A pattern of activity over a set of neurons is presumed to mean something, and there are different ways of coding the same meaning. But different codes have different properties. A code may not be sufficient to uniquely code all the possible things we need to represent. A code could be redundant and have more than one way of representing the same thing. This section continues with our review of the basics of vector and linear algebra by going a little more deeply into the subject. The pay-off will be some mathematics that provides intuition about issues of neural representation. You can think of this as a first lesson in the "psychology of linear algebra".

Basis sets

It is pretty clear that given any vector whatsoever in 8-space, you can specify how much of it gets projected in each of the eight directions specified by the unit vectors v1, v2, ...v8. But you can also build back up an arbitrary vector by adding up all the contributions from each of the component vectors. This is a consequence of vector addition and can be easily seen to be true in 2 dimensions. We can verify it ourselves. Pick an arbitrary vector g, project it onto each of the basis vectors, and then add them back up again:

g = {2,6,1,7,11,4,13, 29} ;

(g.u1) u1  + (g.u2) u2  +(g.u3) u3  +(g.u4) u4  +
(g.u5) u5  +(g.u6) u6  +(g.u7) u7  +(g.u8) u8

Exercise

What happens if you project g onto the normalized Walsh basis set defined by {w1,w2,...}  above, and then add up all 8 components?

(g.w1) w1  + (g.w2) w2  +(g.w3) w3  +(g.w4) w4  +
(g.w5) w5  +(g.w6) w6  +(g.w7) w7  +(g.w8) w8

The projections, g.ui are sometimes called the spectrum of g. This terminology comes from the Fourier basis set used in Fourier analysis. A discrete version of a Fourier basis set is similar to the Walsh set, except that the elements fit a sine wave pattern, and so are not binary-valued.
The orthonormal set of vectors we've defined above  is said to be complete, because any vector in 8-space can be expressed as a linear weighted sum of these basis vectors. The weights are just the projections. If we had only 7 vectors in our set, then we would not be able to express any 8-dimensional vector in terms of this basis set. The seven vector set would be said to be incomplete. A basis set which is orthonormal and complete is very nice from a mathematical point of view. Another bit of terminology is that these seven vectors would not span the 8-dimensional space. But they would span some sub-space, that is of smaller dimension, of the 8-space.
There has been much interest in describing the effective weighting properties of visual neurons in primary visual cortex of higher level mammals (cats, monkeys) in terms of basis vectors. One issue is if the input (e.g. an image) is projected (via a collection of  receptive fields) onto a set of neurons, is information lost? If the set of weights representing the receptive fields of the collection of neurons is complete, then no information is lost.

Linear dependence

What if we had 9 vectors in our basis set used to represent vectors in 8-space? For the u's, it is easy to see that in a sense we have too many, because we could express the 9th in terms of a sum of the others. This set of nine vectors would be said to be linearly dependent. A set of vectors is linearly dependent if one or more of them can be expressed as a linear combination of some of the others.  Sometimes there is an advantage to having an "over-complete" basis set (e.g. more than 8 vectors for 8-space; cf. Simoncelli et al., 1992).
Theorem: A set of mutually orthogonal vectors is linearly independent. 
However, note it is quite possible to have a linearly independent set of vectors which are not orthogonal to each other.  Imagine 3-space and 3 vectors which do not jointly lie on a plane. This set is linearly independent.
If we have a linearly independent set, say of 8 vectors for our 8-space, then no member can be dropped without a loss in the dimensionality of the space spanned.
It is useful to think about the meaning of linear independence in terms of geometry. A set of three linearly independent vectors can completely span 3-space. So any vector in 3-space can be represented as a weighted sum of these 3. If one of the members in our set of three can be expressed in terms of the other two, the set is not linearly independent and the set only spans a 2-dimensional subspace. That is, the set can only represent vectors which lay on a plane in 3-space. This can be easily seen to be true for the set of u's, but is also true for the set of v's.

Thought exercise

Suppose there are three inputs feeding into three neurons in the simple linear network such as defined at the beginning of this lecture. If the weight vectors of the three neurons are not linearly independent, do we lose information?

Basic matrix arithmetic

Definition of a matrix: a list of scalar lists

An mxn matrix has m rows, and n columns. Here is a 3x4 matrix:

MatrixForm[
	Table[W[i,j],{i,1,3},{j,1,4}]]

The matrix can be written in standard form using subscripts as:

And as we have seen, Mathematica represents a matrix as a list of lists:

%//StandardForm

Remember, the symbol % in Mathematica  stands for the most recent output, %% stands for the second to last output, and so forth.

Adding, subtracting and multiplying by a scalar

As with vectors, matrices are added, subtracted, and multiplied by a scalar component by component:

A = {{a,b},{c,d}};
B = {{x,y}, {u,v}};

Let's add A to B:

A+B

Subtracting B from A:

A-B

And if we multiply A by 3:

3 A

Multiplying two matrices

We have already seen how to multiply a vector by a matrix: we replace the ith row of the output vector by the inner product of the ith row of the matrix with the vector.
In order to multiply a matrix A, by another matrix B to get C = AB, we calculate the ijth component of the output matrix by taking the inner product of the ith row of A with the jth column of B:

A.B

Note that AB is not necessarily equal to BA:

B.A

Laws of commutation, association and distribution

Look at the element in the upper left of the matrix BA above--there is no reason, in general, for ax+bu to equal ax + cy. That is, matrix multiplication does not commute.

Apart from commutation for matrix multiplication, the usual laws of commutation, association, and distribution that hold for scalars hold for matrices. Matrix addition and subtraction do commute. Matrix multiplication is associative, so (AB)C = A(BC). The distributive law works too: 
A(B+C) = AB + AC

Non-square matrices

It is not necessary for A and B to be square matrices (i.e. have the same number of rows as columns) to multiply them. But if A is an mxn matrix, then B has to be an nxp matrix in order for AB to make sense. For example, here F is a 3x2 matrix, and G is a 2x4 matrix.

F = {{a,b},{c,d},{e,f}};
G = {{p,q,r,s},{t,u,v,w}};

Dimensions[F]

Dimensions[G]

Because F has 2 columns, and G has 2 rows, it makes sense to multiply G by F:

F.G

However, because the number of columns of G (4)  do not match the number of rows of F (3), G.F is not well-defined:

G.F

Inverse of a Matrix

Dividing a matrix by a matrix: the identity matrix & matrix inverses

The matrix corresponding to 1 or unity is the identity matrix. Like 1, the identity matrix is fundamental enough, that Mathematica  provides a special function to generate n-dimensional identity matrices. Here is a 2x2:

IdentityMatrix[2]

It is easy to show that the identity matrix plays the role for matrix arithmetic that the scalar 1 plays for scalar arithmetic.

How can one divide one matrix, say B, by another, say A? We can divide numbers, x by y, by multiplying x times the inverse of y, i.e. 1/y. So to do the equivalent of dividing B by A, we need to find a matrix Q such that when A is multiplied by Q, we get the matrix equivalent of unity, i.e. the identity matrix. Then "B/A" can be achieved by calculating the matrix product: B.Q.

A = {{a,b},{c,d}};

Mathematica  provides a built-in function to compute matrix inverses:

Q = Inverse[A]

We can test to see whether the product of a A and Q is the identity matrix, but Mathematica  won't go through the work of simplifying the algebra in this case, unless we specifically ask it to.

Q.A

Simplify[Q.A]

Here is a simple numerical example:

B = {{1,-1},{3,2}};
R = Inverse[B]

B.R

Badly conditioned matrices

What if one row is a scaled version of another row? Then the rows are not linearly independent. In this case, the inverse is not defined.

B1 = {{1.5,1},{3,2.0}}

Inverse[B1]

Sometimes the rows are almost, but not quite, linearly dependent (because the elements are represented as approximate floating point approximations to the actual values). Mathematica may warn you that the matrix is badly conditioned. Mathematica  may try to find a solution to the inverse , but you should be suspicious of the solution. In general, one has to be careful of badly conditioned matrices.

Let's try finding the inverseof the following matrix:

B2 = {{-2,-1},{4.00000000000001,2.0}};
Inverse[B2]

B2 . Inverse[B2]

Were we lucky or not?

Determinant of a matrix

There is a scalar function of a matrix called the determinant. If a matrix has an inverse, then its determinant is non-zero. Does B2 have an inverse?

Det[B2]

Why do we get a zero determinant for B2, but not for B1?

Det[B1]

Matrix transpose

We will use the transpose operation quite a bit in this course.  It interchanges the rows and columns of a matrix:

X  =  Table[(Subscript[w, i, j]), {i, 1, 3}, {j, 1, 4}]

You may have noticed that in the Cell menu, you can convert to various cell types.
In StandardForm on the input line, the transpose is written:

Transpose[X];

The output default is TradtionalForm. On the input line, in TraditionalForm, the transpose is written:

(X^T)

What do the ouputs, X and X\(TraditionalForm\`\^T\) look like in StandardForm, MatrixForm?

Getting parts of a matrix

We can pull out the ith row of a matrix in Mathematica by simply writing W[[i]]. For example, the 2nd row of X is:

X[[2]]

What about ith column of a matrix? There is no equally simple way of getting the column of a matrix in Mathematica, but we can use the transpose operation to do it. Transpose[W] [[i]] produces the ith column of matrix W. For example, the 3nd column of X is:

Transpose[X][[3]]

Symmetric matrices

For a square matrix, the diagonal elements remain the same under transpose.
If the transpose of a matrix equals itself, HT = H, H is said to be a symmetric matrix. Symmetric matrices occur quite often in physical systems (e.g. the force on particle i by particle j is equal to the force of j on i). This means that the elements of a symmetric matrix H actually look like they are reflected about the diagonal. We constructed a symmetric matrix in Lecture 5:

H = Table[N[Exp[-Abs[i-j]],1],{i,5},{j,5}]

Neural networks and symmetric connections

Do neural systems have symmetric connections? Real neural networks probably do not in general. Although, when the nature of the processing would not be expected to favor a particular asymmetry, we might expect that there should be symmetric connections on average. We made this assumption when setting up our lateral inhibition weight matrix, as seen above for H. Symmetric matrices have so many nice properties, that neural modelers (especially those from physics backgrounds) find the symmetry assumption almost irresistable. We'll see this later when we study the Hopfield networks. Lack of symmetry can have profound effects on the dynamics of non-linear networks and can produce chaotic trajectories of the state vector.

Outer product of two vectors

We've already seen that the inner product of two vectors produces a scalar. In the next lecture, we will begin our discussion of Hebbian learning in neural networks. In this context, the outer product of an input and output vector will be used to model synaptic modification. Consider two vectors:

f = {f1,f2,f3};
g = {g1,g2,g3};

The outer product is just all of the pairwise products of the elements of f and g arranged in a nice (and special) order:

h = Outer[Times,f,g]

The outer product is also written in traditional row and column format as: f gT

Eigenvectors and eigenvalues

Eigenvectors

An eigenvector, x,  of a matrix, A,  is vector that when you multiply it by A, you get an output vector that points in the same direction as x:
Ax = λx
where λ is a scalar that adjusts the length change of x.
There can't be any more than n distinct eigenvectors for an nxn matrix--and there may be less.

A = {{1,2},{3,4}};

The Mathematica function Eigenvectors[A] returns the eigenvectors of matrix A as the rows of a matrix, which we'll call eig:

eig = Eigenvectors[A]

We can verify that eig[[1]] and A.eig[[1]] lie along the same direction by taking the dot product of the unit vectors pointing in the directions of each:

normalize[x_] := x/Sqrt[x.x];

normalize[eig[[1]]].normalize[A.eig[[1]]];
N[%]

Eigenvalues

The eigenvalues are given by:

Eigenvalues[A]

Eigenvalues and eigenvector elements do not have to be real numbers. They can be complex, that is an element can be the sum of a real and imaginary number. In Mathematica, imaginary numbers are represented by multiples (or fractions) of I, the square root of -1:

Sqrt[-1]
Sqrt[-1]//StandardForm

B = {{1,2},{-3,4}};
Eigenvalues[B]

Linear systems analysis

Introduction

The world of input/output systems can be divided up into linear and non-linear systems. Linear systems are nice because the mathematics that describes them is not only well-known, but also has a mature elegance. On the other hand, it is a fair statement to say that most real-world systems are not linear, and thus hard to analyze...but fascinating if for that reason alone. Scientists were lucky with the limulus eye. That nature is usually non-linear doesn't mean one shouldn't familiarize oneself with the basics of linear system theory. Many times a non-linear system has a sufficiently smooth mapping that it can be approximated by a linear one over restricted ranges of parameter values. 
So exactly what is a "linear system"?
The notion of a "linear system" is a generalization of the input/output properties of a straight line passing through zero. The matrix equation W.x == y is a linear system. This means that if W is  a matrix, x1 and x2 are vectors, and a and b are scalars:
                     W.(a x1 + b x2) = a W.x1 + b W.x2

This is a consequence of the laws of matrix algebra.The idea of a linear system has been generalized beyond matrix algebra.  Imagine we have a box that takes inputs such as f, and outputs g = T[f].

The abstract definition of a linear system is that it satsifies:
	                     T[a f + b g] = a T[f] + b T[g]
where T is the transformation that takes the sum of scaled inputs f, g (which can be functions or vectors) to the sum of the scaled transformation of f and g. The property, that the output of a sum is the sum of the outputs, is sometimes known as the superposition principle for linear systems. The fact that linear systems show superposition is good for doing theory, but as we will see later, it limits the kind of computations that can be done with linear systems, and thus with linear neural network models.

Characterizing a linear system by its response to an orthonormal basis set

Suppose we have an unknown physical system, which we model as a linear system T:

T = Table[Random[],{i,1,8},{j,1,8}];

We would like to make a simple set of measurements that could characterize T in such a way that we could predict the output of T to any input. This is the sort of task that engineers face when wanting to characterize, say a stereo amplifier (as a model linear system), so that the output sound can be predicted for any input sound.What kind of measurements would tell us what T is? Well, we could just "stimulate" the system with cartesian vectors {1,0,0,0,0,0,0,0},{0,1,0,0,0,0,0,0}, and so forth and collect the responses which would be the columns of T. This has two practical problems: 1) for a real physical system, such as your stereo, or a neuron in the limulus eye, this would require stimulating it with a high-intensity audio or light intensity spike, which could damage what you are trying to study; 2) Characterizing the linear system by a matrix T, requires (n^2)numbers, where n is the input signal vector length--and n can be pretty big for both audio and visual systems. Problem 2) has a nice solution when T is symmetric, and even nicer solution if the rows are shifted versions of each other (this is addressed later). Problem 1) can be addressed by showing that we can characterize T with any basis set--so we can pick one that won't blow out the physical system being tested.
The set of Walsh functions we looked at earlier is just one possible set that has the advantage that the elements that contribute to the "energy", i.e. (the square of the length) are distributed across the vector.

Vectorlength[x_] := N[Sqrt[x.x]]

v1 = {1, 1, 1, 1, 1, 1, 1, 1}; w1 = v1/Vectorlength[v1];
v2 = {1,-1,-1, 1, 1,-1,-1, 1}; w2 = v2/Vectorlength[v2];
v3 = {1, 1,-1,-1,-1,-1, 1, 1}; w3 = v3/Vectorlength[v3];
v4 = {1,-1, 1,-1,-1, 1,-1, 1}; w4 = v4/Vectorlength[v4];
v5 = {1, 1, 1, 1,-1,-1,-1,-1}; w5 = v5/Vectorlength[v5];
v6 = {1,-1,-1, 1,-1, 1, 1,-1}; w6 = v6/Vectorlength[v6];
v7 = {1, 1,-1,-1, 1, 1,-1,-1}; w7 = v7/Vectorlength[v7];
v8 = {1,-1, 1,-1, 1,-1, 1,-1}; w8 = v8/Vectorlength[v8];

We have already seen that the set  {(Subscript[w, i])} spans   8-space in such a way that we can easily express any vector as a linear sum of these basis vectors.

g = ∑(g . (Subscript[w, i])) (Subscript[w, i])

So as we saw before, an arbitrary vector, g

g = {2,6,1,7,11,4,13, 29};

is the sum of its own projections onto the basis set:

(g.w1) w1  + (g.w2) w2  +(g.w3) w3  +(g.w4) w4  +
(g.w5) w5  +(g.w6) w6  +(g.w7) w7  +(g.w8) w8

Suppose we now do an "experiment" to find out how T transforms the vectors of our basis set:, and we put all of these transformed basis elements into a new set of vectors newW[[i]]. newW is a matrix for which each row is the response of T to a basis vector.

newW = {T.w1,T.w2,T.w3,T.w4,T.w5,T.w6,T.w7,T.w8};

Note that newW is an 8x8 matrix. So how can we calculate the output of T, given g without actually running the input through T? If we do run the input through T we get:

T.g

But by the principle of linearity, we can also calculate the output by finding the "spectrum" of g as in Problem Set 2, and then scaling each of the transformed basis elements by the spectrum and adding them up:

T . 
        g = T . {∑(g . (Subscript[w, i])) (Subscript[w, i])} = ∑(g . (Subscript[w, i])) 
            T . (Subscript[w, i])

(g.w1) T.w1 + (g.w2) T.w2 + (g.w3) T.w3 + (g.w4) T.w4 + 
(g.w5) T.w5 + (g.w6) T.w6 + (g.w7) T.w7 + (g.w8) T.w8

Of course, we have already done our "experiment", so we know what the transformed basis vectors are, we stored them as rows of the matrix newW. We can calculate what the specturm (g.wi) is, so the output of T is:

(g.w1) newW[[1]] + (g.w2) newW[[2]] + (g.w3) newW[[3]] + 
(g.w4) newW[[4]] + (g.w5) newW[[5]] + (g.w6) newW[[6]] + 
(g.w7) newW[[7]] + (g.w8) newW[[8]]

Same thing in more concise notation

Let  the basis vectors be the rows of a matrix W:

W  =  {w1, w2, w3, w4, w5, w6, w7, w8};

So again, we can project g onto the rows of W, and then reconstitute it in terms of W to get g back again:

(W . g) . W

g . Transpose[W] . newW

The main idea is:  characterize and unknown system T by  its response to orthonormal vectors. As an exercise, you can:

Show that T==Transpose[newW].W,  and that the system output is thus: Transpose[newW].W.g

These new basis vectors do span 8-space, but they are not necessarily orthonormal. Under what conditions would they be orthogonal? What if the matrix T was symmetric, and we deliberately chose the basis set to describe our input to be the eigenvectors of T?

What if the choice of basis set is the set of eigenvectors of T?

We've seen how linearity provides us with a method for characterizing a linear system in terms of the responses of the system to the basis vectors. The problem is that if the input signals are long vectors, say with dimension 40,000, then this set of basis vector responses is really big--1.6x(10^9).
Construct a symmetric matrix transformation, T. Show that if the elements of the basis set are the eigenvectors of T, then the transformation of any arbitrary input vector x is given by:

Where the (Subscript[α, i])are the projections of x onto each eigenvector. Having the eigenvectors of T enables us to express the input and output of T in terms of the same basis set--the eigenvectors. All T does to the input is to scale its projection onto each eigenvector by the eigenvalue for that eigenvector. The set of these eigenvalues is sometimes called the modulation transfer function because it describes how the amplitude of the eigenvectors change as they pass through T.
Linear systems analysis is the foundation of Fourier analysis, and is why it makes sense to characterize your stereo amplifier in terms of frequency response. But your stereo isn't just any linear system--it has the special property that if you input a sound at time t and measure the response, and then you input the same sound again at a later time, you get the same response, except of course that is is shifted in time. It is said to be a shift-invariant system. The eigenvectors of a shift-invariant system are sinusoids. (Recall the eigenvectors of the symmetric matrix in the optional exercises of Lecture 6? They were sinusoids, not just because the matrix was symmetric, but also because each row of the matrix was a shifted version of the previous row--the elements along any given diagonal are identical. This is called a symmetric Toeplitz matrix.)
Sinewave inputs are the eigenvectors of your stereo system. The dimensionality is much higher--if you are interested in frequencies up to 20,000 Hz, your eigenvector for this highest frequency would have least 40,000 elements--not just 8!
This kind of analysis has been applied not only to physical systems, but to a wide range of neural sensory systems. For the visual system alone, linear systems analysis has been used to study the cat retina (Enroth-Cugell and Robson, 1964), the monkey visual cortex, and the human contrast sensivity system as a whole (Campbell and Robson, 1968).
Much empirical analysis has been done using linear systems theory to characterize neural sensory systems, and other neural systems such as those for eye movements. It works wonderfully as long as the linear system approximation holds. And it does do quite well for the lateral eye of the limulus, X-cells and P-cells of the mammalian visual system, over restricted ranges for so-called "simple" cells in the visual cortex, among others. The optics of the simple eye is another example of an approximately linear system. Many non-linear systems can be approximated as linear systems over smooth subdomains.

In summary, if T has n distinct orthogonal eigenvectors, e\(TraditionalForm\`\_i\), and known eigenvalues, (Subscript[λ, i]),  then:
Step 1: Project x onto eigenvectors of T: x.e\(TraditionalForm\`\_i\)
Step 2: Scale each x.e\(TraditionalForm\`\_i\)by the eigenvalue of e\(TraditionalForm\`\_i\): (Subscript[λ, i])x.e\(TraditionalForm\`\_i\)
Step 3: Scale each e\(TraditionalForm\`\_i\) by (Subscript[λ, i])x.e\(TraditionalForm\`\_i\)
Step 4: Sum these up. That's the response of T to x! : (∑_i)(Subscript[λ, i])x . e\(\_i\)

Optional: eigenvectors, matrices and solving algebraic equations

Eigenvectors, eigenvalues: algebraic manipulation

Last time we used the Mathematica function Eigenvectors[] and Eigenvalues[] to produce the eigenvectors and eigenvalues for the matrix equation: Ax = λx. It is worth spending a little time to understand what is being done algebraically. Consider the following pair of equations specified by a 2x2 matrix W acting on xv to produce a scaled version of xv:

W = {{1,2},{2,1}};
xv := {x,y};
lambda xv == W.xv

Finding the eigenvalues is a problem in solving this set of equations. If we eliminate x and y from the pair of equations, we end up with a quadratic equation in lambda:

Eliminate[] & Solve[]

Eliminate[{lambda xv  ==  W . xv, lambda  != 0},  {x, 
        y}]

Solve[-2 lambda + lambda^2 == 3, lambda]

So our eigenvalues are -1 and 3. We can plug these values of lambda into our equations to solve for the eigenvectors:

Solve[{-x == x + 2 y, - y == 2 x + y}, {x,y}]
Solve[{3 x == x + 2 y, 3 y == 2 x + y}, {x,y}]

Reduce

Mathematica is smart enough that we can use Reduce[] to do it all in one line:

Reduce[{lambda xv == W.xv}, {x,y,lambda}]

The eigenvectors are unique only up to a scale factor, so one can choose how to normalize them. For example, we could arbitrarily set x to 1, and then the eigenvectors are: {1,1}, and {1,-1}. Alternatively, we could normalize them to {1/Sqrt[2], 1/Sqrt[2]} and {1/Sqrt[2], -1/Sqrt[2]}.

Eigenvectors[W]

Side note: Solve[] vs. Reduce[]

Solve[] makes assumptions about constraints left unspecified, so the following returns the solution true for any lambda:

Solve[{lambda x == x + 2 y, lambda y == 2 x + y},
	{x,y,lambda}]

Solve[{lambda x == x + 2 y, lambda y == 2 x + y,
lambda !=0, x != 0, y != 0}, {x,y,lambda}]

Reduce[] gives all the possibilities without making specific assumptions about the parameters:

Either of the following forms will work too:

Reduce[lambda xv == W.xv, {xv[[1]],xv[[2]],lambda}]
Reduce[{{lambda x, lambda y} == {x + 2 y, 2 x + y}},
	{x,y,lambda}]

Determinant solution

Ax = λx can be written: (A - Iλ).x, where I is the identity matrix. The interesting values of x that satisfy this equation are the ones that aren't zero. For this to be true, (A - Iλ) must be singular (i.e. no inverse). And this is true if the determinant of (A - Iλ) is zero

References

Simoncelli, E. P., Freeman, W. T., Adelson, E. H., & Heeger, D. J. (1992). Shiftable Multi-scale Transforms. IEEE Trans. Information Theory, 38(2), 587--607.
Strang, G. (1988). Linear Algebra and Its Applications (3rd ed.). Saunders College Publishing Harcourt Brace Jovanovich College Publishers.