MATRICES AND LINEAR ALGEBRA Linear Algebra Matrix manipulation is - - PowerPoint PPT Presentation

MATRICES AND LINEAR ALGEBRA

Linear Algebra

• Matrix manipulation is the original essence of Matlab;

hence the name MATrix LABoratory.

• In this section we will cover the basics of linear algebra, the

ways of using Matlab in this context, and also look into a few applications.

Vectors, Matrices and ND-arrays

>> x = [1,2,3] % vector x = 1 2 3 >> X = [1,2,3; 4,5,6] % matrix X = 1 2 3 4 5 6 >> X3 = cat(3, X, X) % 3D array X3(:,:,1) = 1 2 3 4 5 6 X3(:,:,2) = 1 2 3 4 5 6 >> size(X3) ans = 2 3 2

• We have already seen how

vectors and matrices can be defined; either by tabulating the values or by using functions such as ones, zeros, rand, eye.

• Matlab also supports tables of

higher dimensions; ND-arrays.

• In the example, we define a

3-dimensional table using cat, which catenates the matrices along the specified dimension (here 3).

Matrix Product

>> x = [1, 2, 3, 4]’; % 4-D vector >> T = [1, 0, -1, 0;... 0, 1, 0, 2]; % Projection matrix >> y = T * x y =

• 2

10

−6 −4 −2 2 4 6 8 10 12 −4 −2 2 4 6 8 x1 x2

• Matrix multiplication is the

basic building block of all linear algebra.

• Matrix product is able to

model different kinds of geometric transformations: affine, projective, rotation, shearing, etc.

• Also multitude of algebraic

transformations are most conveniently represented via matrices, e.g., Fourier and Haar.

Example: Fourier Transform

• The Fourier transform is used for detection of frequencies

contained in a signal (e.g., audio).

• The transform is defined as a matrix product Fx, where x is

the signal and the matrix F is         w0

N

w0

N

w0

N

. . . w0

N

w0

N

w−1

N

w−2

N

. . . w−(N−1)

N

w0

N

w−2

N

w−4

N

. . . w−2(N−1)

N

. . . . . . . . . ... . . . w0

N

w−(N−1)

N

w−2(N−1)

N

. . . w−(N−1)2

N

        , with wN = e−2πi/N.

Example: Fourier Transform

• For example, a length N = 4 Fourier transform matrix

becomes:

    e0 e0 e0 e0 e0 e−2πi/4 e−4πi/4 e−6πi/4 e0 e−4πi/4 e−8πi/4 e−12πi/4 e0 e−6πi/4 e−12πi/4 e−18πi/4     =     1 1 1 1 1 i −1 −i 1 −1 1 −1 1 −i −1 i    

• The transform of x = (1, 2, 3, 4) is computed by

    1 1 1 1 1 i −1 −i 1 −1 1 −1 1 −i −1 i         1 2 3 4     =     10 −2 + 2i −2 −2 + 2i    

Fourier Transform in Matlab

% Prepare FT matrix: >> x = cos(2*pi*0.1*(1:N)’); >> powers = (0:N-1)’ * (0:N-1); >> F = exp(-2*pi*i * powers / N); % Transform and plot: >> X = F*x; >> plot(abs(X))

10 20 30 40 50 60 70 80 90 100 −1 −0.5 0.5 1 Time An oscillating signal x 10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 Frequency Fourier transform of x

• The attached code prepares the

Fourier transform matrix and computes the multiplication.

• The changing thing in the matrix

is the exponent, so we first create a matrix of exponents; then pass it through the exp function.

• The frequency of oscillation in

easily measured. The FT in fact produces a shadow component

• n the right, which can be

ignored.

• Note: The command fft is the

proper way to do FT.

Matrix Inversion

>> X = rand(3) X = 0.2348 0.8669 0.3692 0.7350 0.0862 0.6850 0.9706 0.3664 0.5979 >> iX = inv(X) iX =

• 0.9186
• 1.7644

2.5886 1.0383

• 1.0037

0.5089 0.8549 3.4791

• 2.8413

>> X*iX ans = 1.0000

• 0.0000

1.0000

• 0.0000

1.0000

• Matrix inversion attempts to find

the matrix that multiplies the

• riginal matrix to identity.
• "Given matrix X, find Y such that

XY = I."

Eigenvalues and -vectors

6 7 8 9 10 11 12 13 14 15 2 3 4 5 6 7 8 9

• Eigenvalue decomposition is

an advanced matrix algebra topic, that appears in many applications.

• Although the detailed

description is not a topic of this course, the left panel illustrates how eigenvalues are used to characterize high dimensional data.

• The eigenvectors describe the

directions that have the highest variation.

Eigenvalues and -vectors

% Suppose we have data points in X with % size = [200, 2]. % Compute the mean of the point cloud: >> m = mean(X, 1) m = 10.3992 5.7508 % Subtract the mean. We need to extend % the vector m into a 200x2 matrix first % (i.e., copy it 200 times vertically). >> M = repmat(m, [200,1]); % Compute the direction of max. variance % using eig [V,D] = eig(cov(X - M)) V =

• 0.6969

0.7172 0.7172 0.6969 % The max variation is now along direction (x,y) = (0.7172, 0.6969).

• The direction of maximum

variation can be calculated as shown here.

• Note that the mean of the

points has to be subtracted first.

• After that, the axes are solved

by finding the eigenvectors of the covariance matrix of the data.

• This algorithm is called the

principal component analysis (PCA).

Example: Least Squares

• As an example of a typical linear algebra procedure in

Matlab, let us consider least squares fitting of data.

• The Least Squares method is a curve fitting method that
• riginates from 1795, when Gauss invented the method (at

the time he was only 18 years old).

• The method became widely known as Gauss was the only
• ne able to describe the orbit of Ceres, a minor planet in

the asteroid belt between Mars and Jupiter that was discovered in 1801.

Example: Least Squares

5 10 15 20 25 5 10 15 y(n) x(n) LS minimizes the total length of vertical distances

• Least squares minimizes

the distance between the data and the model.

• For example: Find the line

that passes closest to these data points.

• More specifically: Which

line has the smallest average distance to these points.

Example: LS Theory (very briefly)

• Suppose we would like to fit the function f(x) = ax + b

into the data points (x1, y1), (x2, y2), . . . (xN, yN).

• Task is to find best a, b and c.
• This is done by minimizing the sum of distances:

minimize

N

• n=1

(yn − f(xn))2

• In other words

minimize

N

• n=1

(yn − axn − b)2

Example: LS Theory (very briefly)

• In Matlab, we represent the problem with matrices as

follows.

• Minimize the distance between these two vectors:

       ax1 + b ax2 + b ax3 + b . . . axN + b        and        y1 y2 y3 . . . yN       

Example: LS Theory (very briefly)

• The coefficients a and b repeat at each row, so we can write

the task more compactly as follows

• Minimize the distance between these two vectors:

       x1 1 x2 1 x3 1 . . . . . . xN 1       

• X

a b

• c

and        y1 y2 y3 . . . yN       

• y

LS in Matlab

% Prepare matrix X (add the constant term) >> X = [x, ones(size(x))]; % Compute coefficients a and b >> c = inv(X’ * X) * X’ * y c = 0.5372 0.1311 % Compute coefficients a and b % using the shortcut >> c = X \ y c = 0.5372 0.1311

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 Data Best linear fit Best quadratic fit

• It can be shown that the best

coefficients are given by c = (XTX)−1XTy

• This appears so often, that

there is a shortcut in Matlab: c = X \ y.

• In fact, there are also other

shortcuts for this:

• pinv(X) * y
• y’ / X’
• One can also add higher order

terms to fit other polynomials:

X = [x.^3, x.^2, x, x.^0];

Application: Remove Uneven Illumination

• As an application example of the LS,

consider the attached microscope images.

• The illumination is uneven, which

prevents automatic extraction of cells.

• What we do, is to model the illumination

component using a LS model.

• The model attempts to predict the

brightness from the X-Y-coordinates.

• The model can only predict large scale

variations (slowly changing illumination), and not the small details (cells).

• We choose the 2nd order model because

the illumination looks like a paraboloid.

Application: Remove Uneven Illumination

• Target is to create a model f to predict brightness z from

location (x, y); i.e., z ≈ f(x, y).

• The complete model is now

     z1 z2 . . . zN      =      x2

1

y2

1

x1y1 x1 y1 1 x2

2

y2

2

x2y2 x2 y2 1 . . . . . . . . . . . . . . . . . . x2

N

y2

N

xNyN xN yN 1           c1 c2 . . . c6     

• Above, z1, z2, . . . , zN denote the brightness at coordinates

(x1, y1), (x2, y2), . . . , (xN, yN).

• Total 1.3 million points in each vector x, y and z.

How to do it in Matlab

% Load image img = imread(’unevenIllumination.jpg’); img = double(img); % We need all x and y coordinate combinations [y, x] = meshgrid(1:size(img, 1), 1:size(img, 2)); % Vectorize x and y: 1030 x 1300 -> 1339000 x 1 x = x(:); y = y(:); % Try to predict the illumination from X and Y z = img(:); % Prepare the observation matrix X = [x.^2, y.^2, x.*y, x, y, ones(size(x))]; % Calculate coefficients for each component c = X \ z; % Compute the prediction result zHat = X * c; zHat = reshape(zHat, size(img)); % <-- Illumination result = img - zHat; % <-- Details

• Meshgrid is used to

produce all x-y combinations.

• The observation matrix

X contains 6 terms:

• 2nd order terms:

x2, y2, xy.

• 1st order terms: x, y.
• Constant:
• nes(size(x)).
• But you can throw in

whatever components you think will help in prediction: x3, y3, √xy.

Other Linear Algebra Functions in Matlab

det % Determinant svd % Singular value decomposition lu % LU decomposition qr % QR decomposition cholesky % Cholesky decomposition rref % Reduced row echelon form norm % Matrix or vector norm cond % Matrix condition number rank % Matrix rank trace % Sum over the diagonal diag % Extract/create diagonal fliplr % Flip horizontally flipud % Flip vertically transpose % Transpose horzcat % Append matrices left-right vertcat % Append matrices top-down

• The linear algebra

functions in Matlab are too many to list here.

• The list shows some

functions you may have seen in math classes.

• A rule of thumb: Before

you implement, check if the function exists in Matlab or at Mathworks file exchange.1

1http://www.mathworks.com/matlabcentral/fileexchange/

TOOLBOXES

Toolboxes in Matlab

• Toolboxes extend the basic

functionalities of Matlab, and are purchased separately from base Matlab.

• TUT has an extensive license
• f about 50 toolboxes.
• Command ver shows the

installed toolboxes.

• Each toolbox has a general

help of its contents; e.g., help stats lists all functions

• f the Statistics toolbox.
• A list of all helps is given by

help with no arguments.

Signal Processing Toolbox

1000 2000 3000 4000 5000 6000 7000 8000 200 400 600 800 Frequency / Hz Power Signal Before Filtering 1000 2000 3000 4000 5000 6000 7000 8000 0.2 0.4 0.6 0.8 1 Frequency / Hz Attenuation Frequency Response of the Filter Coefficient about zero Coefficient about one 1000 2000 3000 4000 5000 6000 7000 8000 200 400 600 800 Frequency / Hz Power Signal After Filtering Removed Frequencies Removed Frequencies Removed Frequencies

• As time does not allow

in-depth study of many toolboxes, let’s make an excursion to one: The Signal Processing Toolbox.

• Traditional tasks in signal

processing are related to filter design: How to remove certain frequencies from a discrete time signal.

• This can be achieved by

designing a filter, whose frequency response acts as a multiplier for different frequencies.

Convolution

20 40 60 80 100 120 140 160 180 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 n (samples) Amplitude Signal Before Filtering 5 10 15 20 25 30 −0.4 −0.2 0.2 0.4 0.6 0.8 1 n (samples) Amplitude Impulse Response 20 40 60 80 100 120 140 160 180 −0.15 −0.1 −0.05 0.05 0.1 n (samples) Amplitude Signal After Filtering

• The actual filtering operation

is called convolution.

• The filter is defined by an

impulse response h(n).

• Convolution of signal x(n)

with impulse response h(n) is defined as y(n) =

N

• k=0

h(k)x(n − k)

• In other words, the most

recent inputs are weighted by coefficients h(0), . . . , h(N) and summed together.

Filtering in Matlab

% Load example signal. % This imports signal to y and % sample rate to Fs. >> load handel; % Design a filter with max delay 30, and % frequencies above 0.4 * Fs removed >> h = fir1(30, 0.4, ’high’); % Filter: >> z = conv(h, y); % Listen. >> soundsc(z, Fs); % (Command ’sound’ just plays the sound, % but ’soundsc ’ scales to max volume % first).

• Matlab SP toolbox has several

functions for easy design of digital filters, including:

• fir1: Windowed filter design
• fir2: Filter design with

arbitrary response

• firls: Least squares design
• firpm: Equiripple filter design
• All frequencies are relative to

sample rate. In the example, Fs = 8192.

• Thus, the designed filter removes

all energy below 0.4 × 8192 = 3277 Hz.

Image Filtering

• The presented methods can be extended to 2-dimensional

data.

• Below, we compute the 2D FFT of the test image, zero all

low frequencies and apply the inverse 2D FFT.

Image Filtering: Matlab Code

% Load image img = imread(’karpanen.bmp’); % Compute the FFT. The fftshift command % moves the interesting area to the center. X = fftshift(fft2(img)); % Zero a circular region near center % First compute the indices of points near the center. y0 = size(img, 1) / 2; % y coordinate of the center x0 = size(img, 2) / 2; % x coordinate of the center r = 30; % Radius of the area to be zeroed [x, y] = meshgrid(1:size(img,1)); % All xy-coordinate pairs idx = (abs(hypot(x-x0, y-y0)) <= r); % Indices of near-center points X(idx) = 0; % Revert from FFT space to normal space: X = fftshift(X); % Undo the fftshift filtered = real(ifft2(X)); % Take inverse FFT. % The imaginary part may have % very small values, so ignore those.

Other SP Tools

% Plot the spectrogram of an example signal >> [x, Fs] = audioread(’seiska.wav’); >> spectrogram(x, 256, 128, 256, Fs, ’yaxis’);

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1000 2000 3000 4000 5000 6000 7000 8000 Time (s) Frequency (Hz)

• spectrogram: Energy

heatmap on time-frequency axes.

• fft: Fast Fourier Transform
• fft2: Fast Fourier Transform

for images

• dct: Discrete cosine

transform

• decimate: Decrease sample

rate by integer factor

• interp: Increase sample rate