Time Series Analysis With Matlab Tutorials

Tutorial | Time-Series with Matlab

Disclaimer Feel free to use any of the following slides for educational purposes, however kindly acknowledge the source.

Hands-On Time-Series Analysis with Matlab

We would also like to know how you have used these slides, slides, so please send us emails with comments or suggestions. Michalis Vlachos and Spiros Papadimitrio Papadimitriou u IBM T.J. Watson Research Center



About this tutorial tut orial 



The goal of this tutorial is to show you that time-series research (or research in general) can be made fun, when it involves visualizing ideas, that can be achieved with concise programming. programming.

Disclaimer  

We are not affiliated affiliated with Mathworks in any way … but we do like like using Matlab a lot  since it makes our lives easier

Matlab enables us to do that.

Will I be able to use this MATLAB right away after the tutorial?



I am definitely smarter than her , but I am not a time series person, perper-se. I wonder what I gain from this tutorial …


What this tutorial is NOT about Moving averages  Autoregressive Autoregressive Autoregres sive models



Errors and bugs are most likely contained in this tutorial. We might be responsible for some of them.


Overview PART A — The Matlab programming programming environment environment



Forecasting/Prediction  Stationarity 



Seasonality

PART B — Basic mathematics mathematics  Introduction / geometric intuition  Coordinates and transforms  Quantized representations  Non-Euclidean distances PART C — Similarity Search and and Applications  Introduction  Representations  Distance Measures Measures  Lower Bounding  Clustering/Classification/Visualization  Applications



Why does anyone need Matlab? 

Matlab enables the efficient Exploratory Data Analysis (EDA)

“Science progresses through observation” -- Isaac Newton Newton

PART A: Matlab Introduction

Isaac Newton

“The greatest value of a picture is that is forces us to notice what we never expected to see” -- John Tukey

John Tukey


Matlab 

Interpreted Language


History Hist ory of Matlab Matlab (MATr (MATrix ix LABo LABorato ratory) ry) “The most important thing in the programming language is the name. I have recently invented a very good name and now I am looking for a suitable language”. ---- R. Knuth

– Easy code maintenance (code is very compact)

Programmed by Cleve Cleve Moler as an interface for for EISPACK & LINPACK

– Very fast array/vector manipulation – Support for OOP 

Easy plotting and visualization



Easy Integration with other Languages/OS’s



1957: Moler goes to Caltech. Studies numerical Analysis



1961: Goes to Stanford. Works with G. Forsythe on Laplacian Laplacian eigenvalues. eigenvalues.

– Interact with C/C++, COM Objects, DLLs 

– Build in Java support (and compiler) – Multi-Platform Support (Windows, Mac, Linux)



1979: Met with Jack Little in Stanford. Started working on porting it to C



1984: Mathworks Mathworks is founded

Extensive number of Toolboxes

– Image, Statistics, Bioinformatics, etc

1977: First edition of Matlab; 2000 lines of Fortran

– 80 functions (now more than 8000 functions)

– Ability – Ability to make executable files files 

Cleve Moler

com/company/aboutus/founders/origins_of_matlab_wm.html founders/origins_of_matlab_wm.html Video:http://www.mathworks.com/company/aboutus/



Current State of Matlab/Mathw Matlab/Mathworks orks 

Matlab, Simulink, Stateflow



Matlab Matlab version version 7.3, R2006b



Used in variety of industries – Aerospace, – Aerospace, defense, computers, computers, communica communication, tion, biotech



Mathworks still is privately owned



Used in >3,500 Universities, with >500,000 users worldwide



2005 Revenue: >350 M.



2005 Employees: 1,400+



Pricing: – starts from 1900$ (Commercial use), – ~100$ (Student Edition)

Moneyisisbetter betterthan than Money poverty,ififonly onlyfor for poverty, … financialreasons… reasons…… financial



Why does anyone need Matlab? 

Matlab enables the efficient Exploratory Data Analysis (EDA)

“Science progresses through observation” -- Isaac Newton Newton

PART A: Matlab Introduction

Isaac Newton

“The greatest value of a picture is that is forces us to notice what we never expected to see” -- John Tukey

John Tukey


Matlab 

Interpreted Language


History Hist ory of Matlab Matlab (MATr (MATrix ix LABo LABorato ratory) ry) “The most important thing in the programming language is the name. I have recently invented a very good name and now I am looking for a suitable language”. ---- R. Knuth

– Easy code maintenance (code is very compact)

Programmed by Cleve Cleve Moler as an interface for for EISPACK & LINPACK

– Very fast array/vector manipulation – Support for OOP 

Easy plotting and visualization



Easy Integration with other Languages/OS’s



1957: Moler goes to Caltech. Studies numerical Analysis



1961: Goes to Stanford. Works with G. Forsythe on Laplacian Laplacian eigenvalues. eigenvalues.

– Interact with C/C++, COM Objects, DLLs 

– Build in Java support (and compiler) – Multi-Platform Support (Windows, Mac, Linux)



1979: Met with Jack Little in Stanford. Started working on porting it to C



1984: Mathworks Mathworks is founded

Extensive number of Toolboxes

– Image, Statistics, Bioinformatics, etc

1977: First edition of Matlab; 2000 lines of Fortran

– 80 functions (now more than 8000 functions)

– Ability – Ability to make executable files files 

Cleve Moler

com/company/aboutus/founders/origins_of_matlab_wm.html founders/origins_of_matlab_wm.html Video:http://www.mathworks.com/company/aboutus/



Current State of Matlab/Mathw Matlab/Mathworks orks 

Matlab, Simulink, Stateflow



Matlab Matlab version version 7.3, R2006b



Used in variety of industries – Aerospace, – Aerospace, defense, computers, computers, communica communication, tion, biotech



Mathworks still is privately owned



Used in >3,500 Universities, with >500,000 users worldwide



2005 Revenue: >350 M.



2005 Employees: 1,400+



Pricing: – starts from 1900$ (Commercial use), – ~100$ (Student Edition)

Moneyisisbetter betterthan than Money poverty,ififonly onlyfor for poverty, … financialreasons… reasons…… financial



Matlab Mat lab 7.3 

Who needs Matlab? 

R2006b, Released on Sept 1 2006



– Distributed computing

R&D companies for easy application deployment Professors

– Lab assignments

– Better support for large files

– Matlab allows focus on algorithms not on language features

– New optimization Toolbox 

– Matlab builder for Java

Students

– Batch processing of files

• create Java classes from Matlab

• No more incomprehensible incomprehensible perl code!

– Great environment for testing ideas

– Demos, Webinars Webinars in Flash format format

• Quick coding of ideas, then porting to C/Java etc

– Easy visualization

– (http://www.mathworks.com/produ (http://www.mathworks.com/products/matlab/de cts/matlab/demos. mos. html)

– It’s cheap! (for students at least…)


Starting up Matlab

Tutorial | Time-Series with Matlab Personally I'm always ready to learn, although I do not always like being taught.



Dos/Unix like directory navigation



Commands like:

Matlab Matl ab Enviro Environmen nmentt

Sir Winston Churchill

– cd – pwd

Command Window: - type commands - load scripts

– mkdir 

For navigation it is easier to just copy/paste the path from explorer E.g.: cd ‘c:\d ‘c:\docume ocuments\’ nts\’ Workspace: LoadedVariables/Types/Size



Matlab Matl ab Envir Environme onment nt

Matlab Matl ab Enviro Environmen nmentt



Workspace: LoadedVar iables/Types/Size iables/Types/Size

Workspace: LoadedVariables/Types/Size

Help contains a comprehensive comprehensi ve introduction to all functions

Excellent demos and tutorial of the various features and toolboxes



Starting with Matlab

Populating arrays



Everything is arrays



Manipulation of arrays is faster faster than than regular manipulation with for-loops



Plot sinusoid function a = [0:0.3:2*pi] % generate values from 0 to 2pi (with step of 0.3)

b = cos(a) cos(a) % access cos at positions contained in array [a] (y -axis) plot(a,b plot(a,b) ,b) % plot a (x -axis) against bb (y

a = [1 2 3 4 5 6 7 9 10] % define an array

Related: linspace(-100,100,15); % generate 15 values between -100 and 100



Array Access 

2D Arrays

Access array elements >> a(1)

>> a(1:3) ans = 0

ans =



– Let’s define a 2D array 0.3000

0.6000

0 

>> a = [1 2 3; 4 5 6] a =

Set array elements >> a(1) = 100

Can access whole columns or rows

>> a(1,:)

Row-wise access

ans = 1 4

>> a(1:3) = [100 100 100] 100]

2 5

3 6

1

2

>> a(2,2)

>> a(:,1)

ans =

ans =

5

3 Column-wise access

1 4

A good listener is not only popular everywhere, but after a while he gets to know s omething. –Wilson Mizner



Column-wise computation 

Concatenating Concatenatin g arrays

For arrays greater than 1D, all computations happen column-by-column >> a = [1 2 3; 3 2 1] a =



>> max(a) max(a max(a)

2 2

3 1

3

>> mean(a) mean(a mean(a)

>> sort(a)

ans =

ans =

2.0000

2.0000

2.0000

1 3

Row next to row

>> a = [1 2 3]; >> b = [4 5 6]; >> c = [a b]

ans = 1 3

Column-wise or row-wise

2

3

1 3

Column next to column

c = 1

2 2

>> a = [1;2]; >> b = [3;4]; >> c = [a b] c =

2

3

4

5

Row below below row

>> a = [1 2 3]; >> b = [4 5 6]; >> c = [a; b] c =

1 2

6

>> a = [1;2]; >> b = [3;4]; >> c = [a; b] c =

1 4

2 5

3 6

3 4

1 2 3 4

Column below column



Initializing arrays 

Reshaping and Replicating Arrays

Create array of ones [ ones]



>> a = ones(1,3) a = 1

1

>> a = ones(2,2)*5; a = 1

5 5

Changing the array shape [reshape [reshape]]

– (eg, for easier column-wise computation)

5 5

>> a = [1 2 3 4 5 6] ’; % make it into into aa column >> reshape(a,2,3)

>> a = ones(1,3)*inf a = Inf Inf Inf

ans = 1 2



Create array of zeroes [ zeros]



– Good for initializing arrays

0

0

3 4

5 6

Replicating an array [repmat [repmat]] >> a = [1 2 3]; >> repmat(a,1,2)

>> a = zeros(1,4) a = 0

>> a = zeros(3,1) + [1 2 3] ’ a = 1 2 3

0

ans =

1

repmat(X,[M,N]): make [M,N] tiles of X

2

3

>> repmat >> repmat(a,2,1) repmat(a,2,1) ans = 1 2 1 2

3 3

1



Useful Array functions

Useful Array functions



Last element of array [end ]



>> a = [1 3 2 5]; >> a(end)

>> a = [1 3 2 5]; >> a(end -1)

ans =

ans =

2

3

Find a specific element [find ] ** >> a = [1 3 2 5 10 5 2 3]; >> b = find(a==2) b =

2

5 

reshape(X,[M,N]): [M,N] matrix of columnwise version of X

3

Length of array [length]

Length = 4



7

Sorting [ sort] ***

>> length(a) ans =

a=

1

3

2

5

4 

Dimensions of array [size] >> [rows, columns] = size(a) rows = 1 columns = 4

a=

1

3

2

5

s=

1

2

3

5

i=

1

3

2

4

s = columns = 4

1 = s w o r

1

2

3


Visualizing Data and Exporting Figures 

>> a = [1 3 2 5]; >> [s,i]=sort(a)

Use Fisher’s Iris dataset >> load fisheriris

– 4 dimensions, 3 species – Petal length & width, sepal length & width

1

2

3

5

1

3

2

4

i =

5

Indicates the index where the element came from

strcmp, scatter, hold on


Visualizing Data (2D) >> idx_setosa = strcmp(species, ‘setosa setosa’ ’); % rows of setosa data >> idx_virginica = strcmp(species, ‘virginica virginica’ ’); % rows of virginica >> >> setosa = meas(idx_setosa,[1:2]); >> virgin = meas(idx_virginica,[1:2]); >> scatter(setosa(:,1), setosa(:,2)); % plot in blue circles by default >> hold on; red[r [r ] squares squares[s [s] [s] for these >> scatter(virgin(:,1), >> scatter(virgin(:,1), virgin(:,2), ‘rs rs’ ’); % red

– Iris: • virginica/versicolor/setosa meas (150x4 array): Holds 4D measurements ... 'versicolor' 'versicolor' 'versicolor' 'versicolor' 'versicolor' 'virginica' 'virginica' 'virginica' 'virginica‘

idx_setosa ... 1 1 1 0 0 0

... species (150x1 cell array): Holds name of species for the specific measurement

... The world is governed more by appearances rather than realities…--Daniel Webster

An array of zeros and ones indicating the positions where the keyword ‘setosa ‘ setosa ’’was was found

scatter3



Visualizing Data (3D) >> >> >>

Changing Plots Visually

idx_setosa = strcmp(species, ‘setosa’ setosa’); % rows of setosa data idx_virginica = strcmp(species, ‘virginica ’); % rows of virginica idx_versicolor = strcmp(species, ‘versicolor ’); % rows of versicolor

Zoom out

>> setosa = meas(idx_setosa,[1:3]); >> virgin = meas(idx_virginica,[1:3]); >> versi = meas(idx_versicolor,[1:3]); >> scatter3(setosa(:,1), setosa(:,2),setosa(:,3)); % plot in blue circles by default >> hold on; >> scatter3(virgin(:,1), virgin(:,2),virgin(:,3), ‘rs’ red[r ] squares[s squares[s]] for these rs’); % red[r >> scatter3(versi(:,1), >> scatter3(versi(:,1), virgin(:,2),versi(:,3), ‘gx’ gx’); % green x’s

Zoom in

Create line

Computersare are Computers useless.They Theycan can useless. onlygive giveyou you only answers… answers…

Create Arrow 7

6

5

>> grid on; % show grid on axis >> rotate3D on; % rotate with mouse

4

3

Select Object

Add text

2

1 4.5 4

8 3.5

7 3

6 2.5

5 2

7.5

6.5

5.5

4.5 4



Changing Plots Visually

Changing Plots Visually (Example)  Add

titles

 Add

labels on axis



Change color and width of a line

Change tick labels

 Add

grids to axis



Change color of line



Change thickness/ Linestyle



A

Right click C

etc B



Changing Plots Visually (Example)

Changing Figure Properties with Code 

The result …

GUI’s are easy, but sooner or later we realize that coding is faster

>> a = cumsum(randn(365,1)); % random walk of 365 values

If this represents a year’s worth of measurements of an imaginary quantity, we will change:

Other Styles: 3 2 1

• x-axis annotation to months

0 -1

• Axis labels labels

-2 -3

0 3

10

20

30

40

50

60

70

80

90

• Put title in the the figure figure

100

• Include some greek letters in the title just title just for fun

2 1 0 -1 -2 -3

0

10

20

30

40

50

60

70

80

90

100

Real men do it command-line… --Anonymous



Changing Figure Properties with Code

Changing Figure Properties with Code





Axis annotation to months

>> axis tight; % irrelevant but useful... >> xx = [15:30:365]; >> set(gca, ‘xtick’ xtick’,xx)

The result …

The result …





Changing Figure Properties with Code 

Axis annotation to months

Saving Figures Other latex examples:

Axis labels and title

>> set(gca, ’xticklabel ’,[‘ ,[‘Jan’ Jan’; ... ‘Feb’ ‘ Feb’; Mar’ Mar’])

\alpha, \beta, e^{-\alpha} etc



>> title(‘ ‘ \ title( My My measurements ( \epsilon/ pi)’ pi)’)

Matlab allows to save the figures (.fig) for later processing

>> ylabel( ‘Imaginary Quantity ’)

.fig can be later opened through Matlab

>> xlabel( Month ‘ Month of 2005 2005’)


You can always put-off for tomorrow, what you can do today. -Anonymous



Exporting Figures

Exporting figures (code)  You

Export to: emf, eps, jpg, etc



can also achieve the same result with Matlab code

Matlab code: % extract to color eps

print -depsc myImage.eps ; % from commandcommand-line line depsc myImage.eps; print(gcf,’ ’ print(gcf,’-depsc’ depsc’, myImage’ myImage’) % using variable as name



Visualizing Data - 2D Bars

Visualizing Data - 3D Bars data

1

10

2

10 9 8 6 6 3

8

3

6

4

4

colormap

2 0

8 6 6 5 3 2

colormap 7 5 4 4 2 1

64

1

1.0000 1.0000 1.0000 1.0000

2 3 5 6 7

1

bars

3

2

0 0.0198 0.0397 0.0595 0.0794 0.0992

0 0.0124 0.0248 0.0372 0.0496 0.0620 ... 0.7440 0.7564 0.7688 0.7812

0 0.0079 0.0158 0.0237 0.0316 0.0395 0.4738 0.4817 0.4896 0.4975

3

time = [100 120 80 70]; % our data h = bar(time ); % get handle = bar(time); cmap = [1 0 0; 0 1 0; 0 0 1; .5 0 1]; % colors colormap(cmap ); % create colormap

data = [ 10 8 7; 9 6 5; 8 6 4; 6 5 4; 6 3 2; 3 2 1]; bar3([1 bar3([1 2 3 5 6 7], data); c = colormap(gray ); % get colors of colormap c = c(20:55,:); % get some colors colormap(c ); % new colormap

cdata = [1 2 3 4]; % assign colors set(h,'CDataMapping','direct','CData',cdata );



Visualizing Data - Surfaces

Creating .m files

9

1

8

– Script: A series of Matlab commands (no input/output arguments)

2

3 …

10

– Functions: Programs that accept input and return output

1

7

Standard text files



data 10

6

9 10

5 4

1

10

3 2 1 10 8

10 6

8 6

4

4

2

2 0

The value at position x-y of the array indicates the height of the surface

Right click

0

data = [1:10]; data = repmat(data,10,1); % create data surface(data,'FaceColor',[1 1 1], ' Edgecolor ', [0 0 1]); % plot data view(3); grid on; % change viewpoint and put axis lines


cumsum, num2str, save


Creating .m files

Creating .m files The following script will create:



– An array with 10 random walk vectors M editor Double click

– Will save them under text files: 1.dat, …, 10.dat Sample Script

myScript.m

length 100 100 a = cumsum(randn(100,10)); % 10 random walk data of length % number of columns for i=1:size(a,2), data = a(:,i) ; vector of of characters! characters! fname = [num2str(i) ‘.dat’ .dat’]; % a string is aa vector save(fname, ’data’ data’,’ ASCII’ ASCII’); % save each column in a text file end Write this in the M editor…

A random walk time-series

A

cumsum(A)

1

1

2

3

3

6

4

10

5

15

10

5

0

-5

0

10

20

30

40

50

60

70

80

90

100

…and execute by typing the name on the Matlab command line



Functions in .m scripts

Cell Arrays

When we need to:





– Organize our code

– Let’s read the files of a directory

– Frequently change parameters in our scripts

>> f = dir( ‘*.dat’ *.dat’) % read file contents f = 15x1 struct array with fields: name date bytes isdir for i=1:length(f), a{i} = load(f(i).name); N = length(a{i}); plot3([1:N], a{i}(:,1), a{i}(:,2), ... ... ‘r-’, ‘Linewidth ’, 1.5); grid on; 600 pause; 500 cla; 400 end

keyword output argument function name input argument

function dataN = zNorm(data) % ZNORM zNormalization of vector % subtract mean and divide by std

Help Text (help function_name)

if (nargin <1), (nargin <1), % check parameters error(‘ ‘ error( Not Not enough arguments ’); end data = data – mean(data); mean(data); % subtract mean data = data/ std(data); std(data); % divide by std dataN = data;

Cells that hold other Matlab arrays

Function Body

Struct Array 1 2 3 4 5

300 200

ore arguments function [a,b] = myFunc(data, x, y) % pass & return return m more

100 0 1000 1500 500

See also:varargin, varargout

1000 500



Reading/Writing Files

Flow Control/Loops





Load/Save are faster than C style I/O operations

– But fscanf, fprintf can be useful for file formatting or reading non-Matlab files



fid = fopen('fischer.txt', 'wt');





x = 0:.1:1; y = [x; exp(x)]; fid = fopen('exp.txt','w'); fprintf(fid,'%6.2f %12.8f \n',y); fclose(fid); 0.1 1.1052

0.2 1.2214

0.3 1.3499

0.4 0.4 1.4918

0.5 1.64 87

0.6 1.8221

for – Execute statements a fixed number of times

Elements are accessed column-wise (again…)

0 1

while – Execute statements infinite number of times

for i=1:length(species), fprintf(fid, '%6.4f %6.4f %6.4f %6.4f %s \n', meas(i,:), species{i}); end fclose(fid);

Output file:

if (else/elseif) , switch – Check logical conditions



break, continue



return – Return execution to the invoking function

0.7 2.0138

Life is pleasant. Death is peaceful. It’s the transition that’s troublesome. –Isaac Asimov

tic, toc, clear all


For-Loop or vectorization? clear all; tic; for i=1:50000 a(i) a(i) = sin(i); sin(i); end toc

clear all; a = zeros(1,50000); tic; for i=1:50000 a(i) a(i) = sin(i); sin(i); end toc

clear all; tic; i = [1:50000]; a = sin(i); sin(i); toc; toc;



elapsed_time=

– No need for Matlab to resize everytime

5.0070



elapsed_tim e = 0.0200

Time not important…only life important. –The Fifth Element

Functions are faster than scripts


Matlab Profiler 

Find which portions of code take up most of the execution time

– Identify bottlenecks – Vectorize offending code

– Compiled into pseudocode

elapsed_tim e = 0.1400

Pre-allocate arrays that store output results



Load/Save faster than Matlab I/O functions



After v. 6.5 of Matlab there is for-loop vectorization (interpreter)



Vectorizations help, but not so obvious how to achieve many times Time not important…only life important. –The Fifth Element

name date bytes isdir

m e n . a ( ) f 1



Hints &Tips

Debugging 

There is always an easier (and faster) w ay



-- R. Knuth

Not as frequently required as in C/C++

– Set breakpoints, step, step in, check variables values

– Typically there is a specialized function for what you want to achieve 

Beware of bugs in the above code; I have only proved it correct,not tried it

Set breakpoints

Learn vectorization techniques, by ‘peaking’ at the actual Matlab files: – edit [fname], eg – edit mean – edit princomp

Matlab Help contains many vectorization examples



Eitherthis thisman manisis Either deador ormy mywatch watch dead hasstopped. stopped. has


Debugging 


Advanced Features – 3D modeling/Volume Rendering 

Full control over variables and execution path

Very easy volume manipulation and rendering

– F10: step, F11: step in (visit functions, as well) A

B

F10 C



Advanced Features – Making Animations (Example) 

Advanced Features – GUI’s

Create animation by changing the camera viewpoint



Built-in Development Environment

– Buttons, figures, Menus, sliders, etc 3

3

2

2

1 3

0

-1

1

-1

-2

0 0 -1 -2

50



1

0

2

-3

-3 0

0

50

50

2

-3 -1

1 0

1

2

3

4

100

100

-1

0

1

2

3

4

azimuth = [50:100 99: -1:50]; % azimuth range of values for k = 1:length(azimuth), plot3(1:length(a), a(:,1), a(:,2), 'r', 'Linewidth',2); grid on; view(azimuth(k),30); % change change new M(k) M(k) = getframe; getframe; % save the frame end movie(M,20); % play movie movie 20 20 times times See also :movie2avi

Several Examples in Help

– Directory listing

-2

100

0 -1

3

4

– Address book reader – GUI with multiple axis



Advanced Features – Using Java 

Matlab is shipped with Java Virtual Machine (JVM)

Advanced Features – Using Java (Example) 

Stock Quote Query – Connect to Yahoo server



Access Java API (eg I/O or networking)



Import Java classes and construct objects



Pass data between Java objects and Matlab variables

– http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do? objectId=4069&objectType=file disp('Contacting YAHOO server using ...'); disp(['url = java.net.URL (' urlString ')']); end; url = java.net.URL(urlString ); try stream = openStream(url ); ireader = java.io.InputStreamReader(stream ); breader = java.io.BufferedReader(ireader ); connect_query_data = 1; %connect made; catch connect_query_data = -1; %could not connect case; disp(['URL : ' urlString urlString ]); error(['Could not connect connect to to server. server. It It may may be unavailable. Try again later.']); stockdata ={}; return; end



Matlab Toolboxes

In case I get stuck…

 You



help [command] (on the command line) eg. help fft



Menu: help -> matlab help



Matlab webinars



Google groups

can buy many specialized toolboxes from Mathworks

– Image Processing, Statistics, Bio-Informatics, etc

I’vehad hadaawonderful wonderful I’ve evening.But Butthis this evening. wasn’tit… it… wasn’t

– Excellent introduction on various topics 

There are many equivalent free toolboxes too:

– SVM toolbox

– http://www.mathworks.com/company/events/archived_webinars.html?fp

• http://theoval.sys.uea.ac.uk/~gcc/svm/toolbox/

– Wavelets

– comp.soft-sys.matlab

• http://www.math.rutgers.edu/~ojanen/wavekit/

– You can find *anything* here

– Speech Processing

– Someone else had the same

• http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html

problem before you!

– Bayesian Networks • http://www.cs.ubc.ca/~murphyk/Software/BNT/bnt.html



Overview of Part B 1.

Introduction and geometric intuition

2.

Coordinates and transforms

PART B: Mathematical notions Eightpercent percentof of Eight successisisshowing showing success up. up.

3.

4.



Fourier transform (DFT)



Wavelet transform (DWT)



Incremental DWT



Principal components (PCA)



Incremental PCA

Quantized representations 

Piecewise quantized / s ymbolic



Vector quantization (VQ) / K-means

Non-Euclidean distances 

Dynamic time warping (DTW)



Applications

What is a time-series Definition: A Definition: Asequence sequenceof ofmeasurements measurementsover overtime time 

Medicine



Stock Market



Meteorology



64.0 62.8 62.0 66.0 62.0 32.0 86.4 ... 21.6 45.2 43.2 53.0 43.2 42.8 43.2 36.4 16.9 10.0 …

Geology



Astronomy



Chemistry



Biometrics



Robotics

Images

Shapes

Motion capture

ECG

Image

Sunspot

Color Histogram 600

Acer platanoides

400 200 0

5 0

1 0 0

1 0 5

2 0 0

2 0 5

5 0

1 0 0

1 0 5

2 0 0

2 0 5

5 0

1 0 0

1 0 5

2 0 0

2 0 5

400

200

0

800 600

Earthquake

400 200 0

…more to come

Time-Series Salix fragilis

time


Time Series


Time Series

e u l a v

e u l a v

x = (3, 8, 4, 1, 9, 6)

9

x5 x2 x6

8 6

x3 x1 x4

3

4 1

time

time



Sequence of numeric values

– Finite: – N -dimensional vectors/points – Infinite: – Infinite-dimensional vectors



Mean

Variance



Definition:





From now on, we will generally assume zero mean — mean normalization:

Definition:

or, if zero mean, then



From now on, we will generally assume unit variance — variance normalization:



Mean and variance

Why and when to normalize 

Intuitively, the notion of “shape” is generally independent of – Average level (mean) – Magnitude (variance)

variance σ mean μ




Unless otherwise specified, we normalize to zero mean and unit variance


Variance “=” Length

Covariance and correlation



Variance of zero-mean series:



Length of N-dimensional vector (L2-norm):



Definition

or, if zero mean and unit variance, then 

So that:

x2 | | x | |

x1



Correlation and similarity

Correlation “=” Angle



How “strong” is the linear relationship



between x t and y t ? For normalized series,

slope

2.5

residual

Correlation of normalized series:



Cosine law:



So that:

2.5

ρ = -0.23

2 1.5

2

ρ = 0.99

1.5

1

1

0.5

D A 0 C

F E 0 B

-0.5

-0.5

0.5

-1

-1

-1.5

-1.5

-2

-2

-2.5



-2

-1

0

FRF

1

2

-2.5

x

θ x.y -2

-1

0

FRF

1

2

y



Ergodicity

Correlation and distance 

Example

For normalized series,



Assume I eat chicken at the same restaurant every day and



Question: How often is the f ood good? – Answer one:

i.e., correlation and squared Euclidean distance are linearly related.

– Answer two: 

x | | x

| |

θ x.y

Answers are equal

y



Ergodicity

Stationarity

Example

Example



Ergodicity is a common and fundamental assumption, but sometimes can be wrong:



“Total number of murders this year is 5% of the population”



“If I live 100 years, then I wi ll commit about 5 murders, and if I live 60 years, I will commit about 3 murders”



… non-ergodic!



Such ergodicity assumptions on population ensembles is commonly called “racism.”



– Two weeks ago: – Last month: – Last year:




Is the chicken quality consistent? – Last week:



Autocorrelation

ergodic

– “If the chicken is usually good, then my guests today can safely order other things.”

- y

Answers are equal

stationary


Time-domain “coordinates”

Definition:

6 4

3.5 2 1.5





Is well-defined if and only if the series is (weakly) stationary

1

=

-0.5 -2

Depends only on lag ℓ , not time t -0.5

+

4

+ 1.5

+ -2

+ 2

+

6

+ 3.5

+

1



Time-domain “coordinates”

Orthonormal basis 

6 4

Set of N vectors, { e1, e2, …, eN } – Normal: ||ei| | = 1, for all 1 ≤ i ≤ N

3.5

– Orthogonal: ei¢e j = 0, for i ≠ j

2 1.5

1

=

-0.5



-2

Describe a Cartesian coordinate system – Preserve length (aka. “Parseval theorem”) – Preserve angles (inner-product, correlations)

x1 -0.5

£ e1

+ x42

£ e2

x3 + 1.5

£ e3

+ x-24

£ e4

+ x25

£ e5

+ x66

£ e6

x7 + 3.5

£ e7

+ x18

£ e8



Orthonormal basis

Orthonormal bases

Note that the coefficients x i w.r.t. the basis { e1, …, eN } are the corresponding “similarities” of x to each basis vector/series:





– Each coefficient is simply the value at one time instant



4

– Time/scale (wavelets)

3.5 2 1

-0.5

What other bases may be of interest? Coefficients may correspond to: – Frequency (Fourier)

6

1.5

The time-domain basis is a trivial tautology:

-0.5

=

+

-2

4

e1

x

– Features extracted from series collection (PCA)

+ … e2

x2



Frequency domain “coordinates”

Time series geometry

Preview

Summary



6

– Series / vector

4

3.5

– Mean: “average level”

2 1.5

– Variance: “magnitude/length”

1

=

-0.5

– Correlation: “similarity”, “distance”, “angle” – Basis: “Cartesian coordinate system”

-2

5.6

- 4.9

Basic concepts:

+ -2.2

+

0

+ 2.8

+ -3

+

0

+ 0.05



Time series geometry

Overview

Preview — Applications

 

The quest for the right basis… Compression / pattern extraction

1.


2.

Coordinates and transforms 





– Similarity / distance



Incremental DWT

– Indexing




– Clustering



Incremental PCA

– Filtering

3.

– Forecasting

Quantized representations

– Periodicity estimation – Correlation

4.











Frequency

Frequency and time

. = 0 

One cycle every 20 time units (period)



= 8 period 20? Why period is the



It’s not 8, because its “similarity” (projection) to a period-8 series (of the same length) is zero.



Frequency and time

Frequency and time

.

. = 0

= 0

period = 10

period = 40



Why is the cycle 20?



Why is the cycle 20?







…and so on



Frequency

Frequency

Fourier transform - Intuition

Fourier transform - Intuition



To find the period, we compared the time series with sinusoids of many different periods





Therefore, a good “description” (or basis) would consist of all these sinusoids



This is precisely the idea behind the discrete Fourier transform

Technical details: – We have to ensure we get an orthonormal basis – Real form: sines and cosines at N /2 different frequencies – Complex form: exponentials at N different frequencies

– The coefficients capture the similarity (in terms of amplitude and phase) of the series with sinusoids of different periods



Fourier transform

Fourier transform

Real form

Real form — Amplitude and phase



For odd-length series,



Observe that, for any f k , we can write

where 

The pair of bases at frequency f k are

are the amplitude and phase, respectively. plus the zero-frequency (mean) component



Fourier transform

Fourier transform

Real form — Amplitude and phase

Complex form





The equations become easier to handle if we allow the series and the Fourier coefficients X k to take complex values:



Matlab note: fft omits the scaling factor and is not unitary—however, ifft includes an scaling factor, so always ifft(fft(x)) == x .

It is often easier to think in t erms of amplitude r k and phase k – e.g.,

1

0.5

0

5

-0.5

-1 0

10

20

30

40

50

60

70

80



Fourier transform

Other frequency-based transforms

Example



2 P 1 B G 0

– Matlab: dct / idct

1 frequency

-1 2 P 1 B G 0

2 frequencies

Discrete Cosine Transform (DCT)



Modified Discrete Cosine Transform ( MDCT)

-1 2 P 1 B G 0

3 frequencies

-1

2 P 1 B G 0

5 frequencies

-1

2 P 1 B G 0

10 frequencies

-1

2 P 1 B G 0

20 frequencies

-1


Overview 1.


2.


3.









Incremental DWT






Incremental PCA


Frequency and time

e.g., .

4.

Piecewise quantized / symbolic







Frequency and time 

Fourier is successful for summarization of series with a few, stable periodic components



However, content is “smeared” across frequencies when there are – Frequency shifts or jumps, e.g.,

– Discontinuities (jumps) in time, e.g.,

≠ 0

period = 10

≠ 0

.


period = 20

etc…



What is the cycle now?



No single cycle, because the series isn’t exactly similar with any series of the same length.


Frequency and time 

If there are discontinuities in time/frequency or frequency shifts, then we should seek an alternate “description” or basis



Main idea: Localize bases in time – Short-time Fourier transform (STFT) – Discrete wavelet transform (DWT)



Frequency and time

Frequency and time

Intuition

Intuition





What if we examined, e.g., eight values at a time?



Can only compare with periods up to eight.

What if we examined, e.g., eight values at a time?

– Results may be different for each group (window)



Frequency and time

Wavelets

Intuition

Intuition



Main idea – Use small windows for small periods • Remove high-frequency component, then

– Use larger windows for larger periods • Twice as large

– Repeat recursively



Can “adapt” to localized phenomena



Fixed window: short-window Fourier (STFT)





Technical details – Need to ensure we get an orthonormal basis

– How to choose window size?

Variable windows: wavelets



Wavelets

Wavelets

Intuition

Intuition — Tiling time and frequency

y c n e u q e r F

Time

) y c n e u q e r f ( e l a c S

y c n e u q e r F

Time

Fourier, DCT, …

) y c n e u q e r f ( e l a c S

y c n e u q e r F

Time

Time

STFT

Wavelets



Wavelet transform

Wavelet transform

Pyramid algorithm

Pyramid algorithm

High pass

High pass Low pass Low pass



Wavelet transform

Wavelet transform

Pyramid algorithm

Pyramid algorithm

w1

High pass

x ≡ w0

High pass

Low pass

v1

Low pass

Low pass


Wavelet transforms

General form

Other filters —examples

A high-pass / low-pass filter pair

v2

High pass Low pass

w3 v3


Wavelet transforms 

w2

High pass

Haar (Daubechies-1)

– Example: pairwise difference / average (Haar) – In general: Quadrature Mirror Filter (QMF) pair • Orthogonal spans, which cover the entire space

– Additional requirements to ensure orthonormality of overall transform… 

Daubechies-2

Daubechies-3

Use to recursively analyze into top / bottom half of frequency band Daubechies-4

Wavelet filter , or Mother filter (high-pass)

Scaling filter , or Father filter (low-pass)

B W e o t t e r s r e f r e t i m q e u e l n o c c a y i l i z s a o l a t i t o i n o n



Wavelets

Wavelets

Example

Example

Wavelet coefficients (GBP, Haar)

500

1000

1500

2000

2500

500

1

1 0 200

1

400

600

800

1000

1200

200

300

400

500

600

100

150

200

250

300

600

800

1000

200

300

400

500

40

60

80

100

120

140

160

50

100

150

200

250

300

20

5

40

60

80

100

120

140

20

30

40

50

60

70

80

5

10

20

30

40

50

60

70

4 D

0.4 0.2 0 -0.2 -0.4

-5 5

10

15

20

25

30

35

40

20

5 D

0.5 0 -0.5

-20 10

15

20

25

30

35

40

5

10

15

20

25

30

35

40

45

5

10

15

20

25

30

35

40

2500

500

1000

1500

2000

2500

1000

2000

2500

500

1000

1500

2000

2500

500

1000

1500

2000

2500

500

1000

1500

2000

2500

500

1000

1500

2000

2500

500

1000

1500

2000

2500

500

1000

1500

2000

2500

500

1000

1500

2000

2500

0.2 0 -0.2 -0.4 -0.6

1500

2000

2500

2000

2500

0.2 0 -0.2 -0.4 0.5

500

1000

1500

2000

2500

500

1000

1500

2000

2500

0.5 0 -0.5

2 1 6 A 0 -1

0.5 0 -0.5 2 1 0 -1

45

500

1000

1500

2000

2500


Multi-resolution analysis (GBP, Daubechies-3)



Wavelet GUI: wavemenu



Single level: dwt / idwt



Multiple level: wavedec / waverec

Analysis levels are orthogonal, 1000

1500

2000

2500

500

1000

1500

2000

2500

500

1000

1500

2000

2500

500

1000

1500

2000

2500

1000

2500 1500

2000

2500

1000

1500

2000

2500

1000 2000

1500

2000

2500

500

1000

1500

2000

2500

500

1000

1500

2000

2500

Di¢D j = 0, for i ≠ j 0 -0.2 -0.4 -0.6

500 0.2 0 -0.2

1000 2

1500

2000

2500

0.2 Haar analysis: simple, piecewise constant 0 -0.2 -0.4 -0.6

1 500

0.4 0.2 3 D -0.20 -0.4

0

1000

1500

2000

2500

-1

500

1000

500 1500

1000 2000

2500

0.2 0 -0.2 -0.41500

2000 500

– wmaxlev 

Wavelet bases: wavefun

0.2

0 Daubechies-3 analysis: less artifacting

2 500

0.5 5 0 D -0.5

1 1000

1500

2000

2500

0

-0.2 -0.4 0.5

500

0

-1

-0.5

500

1000

1500

500

1000

1500

500

2000

2500 1000

2000

2500

0.5 0 -0.5

0.5 0 -0.5

1500 500

2500

2 1 0 -1 500

1000

1500

2000

2500



Other wavelets

More on wavelets



Only scratching the surface…





Wavelet packets – All possible tilings (binary) – Best-basis transform



1500

2 1 0 -1

500

2 1 6 0 A -1

1000

0

Matlab

0.4 0.2 0 -0.2 -0.4

1500

1000

Example

0.1 0 1 -0.1 D -0.2 -0.3

6 D

2000

Wavelets

Multi-resolution analysis (GBP, Haar)

4 D

1500

Wavelets

P 2 B 1 G 0 -1

500

-0.5


2 D

1000

500

6 D

0 5

500

500

80

0

0

2500

0.2 0 -0.2 -0.4

160

0 10

2000

0.4 0.2 3 D -0.20 -0.4

-2

-5

1500

0.2 0 -0.2

600

0 20

1000

0 -0.2 -0.4 -0.6

1200

2

-10

-20

100

1

0

20

400

0 50

500

2500

2 D

-1

5 0 W -5

6 V

200

-1 100

2

10

2000

0

2

6 W

1500

2 1 0 -1

0.1 0 1 -0.1 D -0.2 -0.3

1

3 0 W -2

5

1000

-1

-1

2 0 W -1

Multi-resolution analysis (GBP, Daubechies-3)

P 2 B 1 G 0 -1

2 1 0 -1

1 0 W

4 0 W -2

Multi-resolution analysis (GBP, Haar)

Wavelet coefficients (GBP, Daubechies-3)

P 2 B 1 G 0 -1

Overcomplete wavelet transform (ODWT), aka. maximum-overlap wavelets (MODWT), aka. shiftinvariant wavelets

Signal representation and compressibility Partial energy (GBP)

100

90

80

80

70

70

) y g 60 r e n e 50 % ( y t i 40 l a u Q 30

) y g 60 r e n e 50 % ( y 40 t i l a u Q 30 Time FFT Haar DB3

20

Further reading: 1. Donald B. Percival, Andrew T. Walden, Wavelet Methods for Time Series Analysis , Cambridge Univ. Press, 2006. 2. Gilbert Strang, T ruong Nguyen, Wavelets and Filter Banks , Wellesley College, 1996. 3. Tao Li, Qi Li, Shenghuo Zhu, Mitsunori Ogihara, A Survey of Wavelet Applications in Data Mining , SIGKDD Explorations, 4(2), 2002.

Partial energy (Light)

100

90

10

Time FFT Haar DB3

20

10

0

0 0

2

4

6

Compression(% coefficients)

8

10

0

5

10

Compression (% coefficients)

15



More wavelets

Overview



1.


2.




Keeping the highest coefficients minimizes total error (L2-distance) Other coefficient selection/thresholding schemes for different error metrics (e.g., maximum per-instant error, or L -dist.) – Typically use Haar bases

3.

Further reading: 1. Minos Garofalakis, Amit Kumar, Wavelet Synopses for General Error Metrics , ACM TODS, 30(4), 2005. 2.Panagiotis Karras, Nikos Mamoulis, One-pass Wavelet Synopses for Maximum-Error Metrics, VLDB 2005.


4.









Incremental DWT






Incremental PCA









Wavelets

Wavelets

Incremental estimation




Wavelets

Wavelets





Wavelets

Wavelets



post-order traversal



Wavelets

Overview




Forward transform

:

– Post-order traversal of wavelet coefficient tree

1.


2.


– O(1) time (amortized) – O(logN) buffer space (total) 

Inverse transform:

constant factor: filter length

– Pre-order traversal of wavelet coefficient tree – Same complexity

3.

4.









Incremental DWT






Incremental PCA








Time series collections



Overview

Time series collections



Fourier and wavelets are the most prevalent and successful “descriptions” of time series.





Next, we will consider collections of M time series, each of length N .

Some notation:

– What is the series that is “most similar” to all series in the collection? – What is the second “most similar”, and so on…

values at time t , xt i-th series, x(i)



Principal Component Analysis

Principal component analysis

Example

Exchange rates (vs. USD) D 2 U 0 A -2

1

= 48%

0

-0.05 0.05

F 2 E B 0 -2

2 u U

D 2 A 0 C -2

U u

2

0.05

3

3

500

M 2 E 0 D -2

1000

1500

2000

2500

P 2 S 0 E -2

50

SEK

+ 11% = 92%

0

-0.05 0.05 4 0 4 U u -0.05

First two principal components

2 D A 0 C -2

+ 33% = 81%

0

-0.05

F 2 R 0 F -2

2 Y P 0 J -2

(μ ≠ 0)

Principal components 1-4

0.05

1

U u

40 P 2 B 0 G -2

+ 4% = 96%

AUD

30

F 2 R F 0 -2

Time

x(2) = 49.1u1 + 8.1u2 + 7.8u3 + 3.6u4 +

“Best” basis : { u1, u2, u3, u4 } 1

2 G L 0 N -2 L 2 Z 0 N -2

20 2 , i

υ

Coefficients of each time series w.r.t. basis { u1, u2, u3, u4 } :

F 2 E 0 B -2

10

NZL CHF

0

P 2 S E 0 -2

2 G L 0 N -2

-10

K 2 E S 0 -2

M 2 E 0 D -2

-20

2 Y P J 0 -2

F 2 H C 0 -2 P 2 B 0 G -2

-30 500

1000

Time

1500

2000

-20

-10

0

10

20

30

40


60




Matrix notation — Singular Value Decomposition (SVD)


X = U VT X

50

υ i,1

2500

X = U VT

U

X

U

VT

VT v’1

x(1) x(2)

x( M )

=

u1 u2

uk

.

1

2

3

x(1) x(2)

M

x( M )

=

u 1 u2

uk

.

v’2 1

2

3

N

v’k

time series

coefficients w.r.t. basis in U (columns)

basis for time series

time series


basis for measurements (rows) coefficients w.r.t. basis in U (columns)




Principal component analysis


Properties — Singular Value Decomposition (SVD)

X = U VT X



V are the eigenvectors of the covariance matrix XTX, since



U are the eigenvectors of the Gram (inner-product) matrix XXT, since

U VT σ 1

x(1) x(2)

x( M )

=

u1 u2

uk

.

v1

σ 2

. σ k

scaling factors time series


v2

vk

basis for measurements (rows) Further reading: 1. Ian T. Jolliffe, Principal Component Analysis (2nd ed), Springer, 2002. 2. Gilbert Strang, Linear Algebra and Its Applications (4th ed), Brooks Cole, 2005.



Kernels and KPCA 

Multidimensional scaling (MDS)

What are kernels?



Exchange rates

– Usual definition of inner product w.r.t. vector coordinates is x¢y = ∑i x y i i

SEK CAD AUD

– However, other definitions that preserve

FRF DEMBEF NLG CHF

NZL



the fundamental properties are possible 

Why kernels?

Kernels are still “Euclidean” in some sense – We still have a Hilbert (inner-product) space, even though it may not be the space of the original data

ESP GBP

JPY

For arbitrary similarities, we can still find the eigendecomposition of the similarity matrix – Multidimensional scaling (MDS)

– We no longer have explicit “coordinates”

– Maps arbitrary metric data into a

• Objects do not even need to be numeric

low-dimensional space

– But we can still talk about distances and angles

Exchange rates SEK

– Many algorithms rely just on these two concepts

ESP GBP

CAD AUD

NZL

Further reading: 1. Bernhard Schölkopf, Alexander J. Smola, Learning with Kernels: Support Vector Machines, Regularization, Optimization and Beyond , MIT Press, 2001.

JPY


Principal components

FRF DEMBEF NLG CHF


Matlab

PCA on sliding windows



pcacov





princomp

Empirical orthogonal functions (EOF), aka. Singular Spectrum Analysis (SSA)



[U, S, V] = svd(X)



If the series is stationary, then it can be shown that



[U, S, V] = svds(X, k)

– The eigenvectors of its autocovariance matrix are the Fourier bases – The principal components are the Fourier coefficients

Further reading: 1. M. Ghil, et al., Advanced Spectral Methods for Climatic Time Series, Rev. Geophys., 40(1), 2002.


Overview 1.


2.


3.

4.









Incremental DWT






Incremental PCA









Principal components Incremental estimation



PCA via SVD on X – Singular values

N M

— recap:

2  k £k (diagonal)

• Energy / reconstruction accuracy

– Left singular vectors U 2  N £k • Basis for time series • Eigenvectors of Gram matrix XXT

– Right singular vectors V 2  M £k • Basis for measurements’ space • Eigenvectors of covariance matrix XTX






Incremental estimation — Example



N M

PCA via SVD on X

— recap: First series

k £k

– Singular X values

2 U (diagonal) VT

• Energy / reconstruction accuracy

v1

σ 1

£k – Left singular vectors U 2  N . x(1) x(2)

x( M )

=

u1 u2

30oC

.

σ 2

uk

• Basis for time series • Eigenvectors of Gram matrix XXT

v2 ) 1 (

x

σ k

vk

20oC

– Right singular vectors V 2  M £k

s e i r e S

• Basis for measurements’ space • Eigenvectors of covariance matrix XTX

First three values Other values







Correlations:

First series 30oC

Second series ) 2 (

Let’s take a closer look at the first three measurementpairs…

) 2 (

x

20oC

30oC

x

s e i r e S

s e i r e S 20oC

20oC






) 2 (

x

s e i r e S 20oC




30oC

30oC

Series x(1)

First three lie (almost) on a line in the space of t measurement-pairs… n n e p o m c o a l i p c i n p r Î O( M ) numbers for = e t the slope, and s f o f Î One number for

Other pairs also follow the same pattern: they lie (approximately) on this line

30oC

) 2 (

x

s e i r e S 20oC

each measurementpair (offset on line = PC) 20oC

Series x(1)

30oC


20oC

Series x(1)

30oC







Incremental estimation —Example (update)

For each new point 30oC

error



Project onto current line



Estimate error


) 2 (

) 2 (

s e i r e S 20oC

s e i r e S 20oC

x

error

x

20oC

30oC

20oC

New value

Series x(1)




Estimate error



Rotate line in the direction of the error and in proportion to its magnitude

Î

O (M ) time

30oC

Series x(1)



Incremental estimation —Example (update)



) 2 (

x

s e i r e S 20oC






Estimate error



Rotate line in the direction of the error and in proportion to its magnitude



The “line” is the first principal component (PC) direction



This line is optimal: it minimizes the sum of squared projection errors

30oC

Series x(1)





Incremental estimation —Update equations

Incremental estimation — Complexity

O(M k ) space (total) and time (per tuple), i.e.,

For each new point xt and for j = 1, …, k :  

y j := 2

j

v Tx j

(proj. onto v j )

t

j

+ y j 2



e j := x – y j w j



v j



xt

v j + (1/

xt – y j v j

(energy

j -th eigenval.)

(error) 2)

j

y j e j

New value



20oC



(update estimate) (repeat with remainder)

xt

v1 updated

e 1

v1

y1



Independent of # points



Linear w.r.t. # streams ( M )



Linear w.r.t. # principal components (k )




Overview

Incremental estimation — Applications



Incremental PCs (measurement space)

1.


– Incremental tracking of correlations

2.


– Forecasting / imputation – Change detection

3.

Further reading: 1. Sudipto Guha, Dimitrios Gunopulos, Nick Koudas, Correlating synchronous and asynchronous data streams , KDD 2003. 2. Spiros Papadimitriou, Jimeng Sun, Christos Faloutsos, Streaming Pattern Discovery in Multiple Time-Series , VLDB 2005. 3. Matthew Brand, Fast Online SVD Revisions for Lightweight Recommender Systems , SDM 2003.

4.









Incremental DWT






Incremental PCA









Piecewise constant

Piecewise constant (APCA) 

Example

So far our “windows” were pre-determined

APCA (k=10)

– DFT: Entire series

2 1

– STFT: Single, fixed window

0

– DWT: Geometric progression of windows 


-1

500

Within each window we sought fairly complex patterns (sinusoids, wavelets, etc.)

1000

1500

2000

2500

2000

2500

2000

2500

2000

2500

2000

2500

2000

2500

APCA (k=21) 2 1 0



Next, we will allow any window size, but constrain the “pattern” within each window to the simplest possible (mean)

-1

500

1000

1500

APCA (k=41) 2 1 0 -1

500


1500


Piecewise constant

Piecewise constant (APCA) 

1000

Example

Divide series into k segments with endpoints

APCA (k=10) 2

– Constant length: PAA

1 0

– Variable length: APCA 

-1

Single-level Haar smooths,

if t j+1-t j = 2ℓ , fortheir all 1 ≤ j ≤ k Represent all points within one segment with 1 m j k average j , , thus minimizing 2 1 0 -1

500

1000

1500

2000

500

1000

1500

APCA (k=21) / Haar(level 7, 21 coeffs) 2 1 0 -1

500

1000

1500

APCA (k=41) / Haar (level 6, 41 coeffs) 2

Further reading: 1. Kaushik Chakrabarti, Eamonn Keogh, Sharad Mehrotra, Michael Pazzani, Locally Adaptive Dimensionality Reduction for Indexing Large Time Series Databases , TODS, 27(2), 2002.

1 0 -1

500

1000

1500



Piecewise constant

Piecewise constant

Example

Example

APCA (k=10)

APCA (k=10)

2

2

1

1

0

0

-1

-1

500

1000

1500

2000

2500

500

APCA (k=21) / Haar (level 7, 21 coeffs) 2

2

1

1

0

0

-1

-1

500

1000

1500

2000

2500

500

1500

2000

2500

1000

1500

2000

2500

APCA (k=15) / Daubechies-3 (level 7, 15 coeffs)

APCA / Haar(top-21 out of 7 levels) 2

2

1

1

0

0

-1

-1 500

1000

1500

2000

500

2500


k/h-segmentation 

Again, divide the series into k segments (variable length)



For each segment choose one of h quantization levels to represent all points – Now, m j can take only h ≤ k possible values



1000

APCA (k=21) / Haar(level 7, 21 coeffs)

1000

1500

2000

2500


Symbolic aggregate approximation (SAX) 

Quantization of values



Segmentation of time based on these quantization levels



More in next part…

APCA = k /k -segmentation (h = k )

Further reading: 1. AristidesGionis, Heikki Mannila, Finding Recurrent Sources in Sequences, Recomb 2003.



Overview

K-means / Vector quantization (VQ)

1.




2.


3.

4.









Incremental DWT






Incremental PCA








APCA considers one time series and – Groups time instants – Approximates them via their (scalar) mean



Vector Quantization / K-means applies to a collection of M time series (of length N ) – Groups time series – Approximates them via their (vector) mean



K-means

K-means 

Partitions the time series x(1), …, x(M ) into groups, I j , for 1 j k .



All time series in the j -th group, 1 j represented by their centroid, m j .



Objective is to choose m j so as to minimize the overall squared distortion,

m2

k

k , are

m1 1-D on values + contiguity requirement: APCA



K-means

K-means

m2



Objective implies that, given I j , for 1

j

k ,

i.e., m j is the vector mean of all points in cluster j . m1


K-means


K-means Example

1. Start with arbitrary cluster assignment. 2. Compute centroids.

Exchange rates PCs

50

ESP SEK

3. Re-assign to clusters based on new centroids.

40

4. Repeat from (2), until no improvement.

30

0.05 0

GBP

-0.05 0.05

CAD

0 -0.05

AUD

k = 2 1 0 -1

20 2 , i

Finds local optimum of D .

υ 10

DEM NZL

σ ≠ 1

k = 4

CHF

0

Matlab: [idx, M] = kmeans(X’, k)

2 1 0 -1

FRF BEF NLG

2 1 0 -1 -10

2 1 0 -1 2 0 -2

-20

JPY

-30

-20

-10

0

10

20

υ i,1

30

40

2 0 -2 50

60

σ ≠ 1



K-means in other coordinates

K-means and PCA



An orthonormal transform (e.g., DFT, DWT, PCA) preserves distances.



K-means can be applied in any of these “coordinate systems.”



Can transform data to speed up distance computations (if N large)

Further reading: 1. Hongyuan Zha, Xiaofeng He, Chris H.Q. Ding, Ming Gu, Horst D. Simon,Spectral Relaxation for K-means Clustering , NIPS 2001.



Overview


1.




2.


So far we have been discussing shapes via various kinds of “features” or “descriptions” (bases)



However, the series were always fixed



Dynamic time warping:

3.

4.









Incremental DWT






Incremental PCA

– Allows local deformations (stretch/shrink)






– Can thus also handle series of different lengths





Dynamic time warping (DTW) 


Euclidean (L2) distance is



Each cell c = (i , j ) is a pair of indices whose corresponding values will be compared, (x i – y j )2, and included in the sum for the distance



Euclidean path:

or, recursively,



Dynamic time warping distance is ] j : 1 [ y

– i = j always – Ignores off-diagonal cells

where x1:i is the subsequence (x 1, …, x i )

shrink x / stretch y stretch x / shrink y

x[1:i]



Dynamic time-warping


Fast estimation



DTW allows any path



Examine all paths: shrink x / stretch y (i-1, j)

] j : 1 [ y

(i, j) (i, j)



Standard dynamic programming: O (N 2)



Envelope-based technique – Introduced by [Keogh 2000 & 2002]

s t r e t c h

– Multi-scale, wavelet-like technique and formalism by [Salvador et al. 2004] and, independently, by [Sakurai et al. 2005]

x

(i-1, j- 1) (i, j-1) /

s h r i n k



x[1:i]

y

Standard dynamic programming to fill in table—top-right cell contains final result

Further reading: 1. Eamonn J. Keogh, Exact Indexing of Dynamic Time Warping , VLDB 2002. 2. Stan Salvador, Philip Chan, FastDTW: Toward Accurate Dynamic Time Warping in Linear Time and Space , TDM 2004. 3. Yasushi Sakurai, Masatoshi Yoshikawa, Christos Faloutsos, FTW: Fast Similarity Under the Time Warping Distance , PODS 2005.



Dynamic time warping

Non-Euclidean metrics

Fast estimation —Summary



Create lower-bounding distance on coarser granularity, either at



More in part 3

– Single scale – Multiple scales ] j : 1 [ y



Use to prune search space

x[1:i]



Timeline of part C – Introduction – TimeTime-Series Representations

PART C: Similarity Search and Applications

– Distance Measures – Lower Bounding – Clustering/Classification/Visualization – Applications



Applications (Image Matching)

Applications (Shapes) Cluster 1 Recognize type of leaf based on its shape

Many types of data can be converted to time-series Image

Ulmus carpinifolia

Color Histogram 200 0

5 0

1 0 0

1 0 5

2 0 0

2 0 5

5 0

1 0 0

1 0 5

2 0 0

2 0 5

5 0

1 0 0

1 0 5

2 0 0

2 0 5

Salix fragilis

Tilia

Quercus robur

Convert perimeter into a sequence of values

600 400

Acer platanoides

Cluster 2

400

200

0

800 600 400 200 0

Time-Series Special thanks to A. Ratanamahatana & E. Keogh for the leaf video.


Applications (Motion Capture) Motion-Capture (MOCAP) Data (Movies, Games)

– Track position of several joints over time – 3*17 joints = 51 parameters per frame


Applications (Video) Video-tracking / Surveillance

– Visual tracking of body features (2D time-series) – Sign Language recognition (3D time-series) Video Tracking of body feature over time (Athens1, Athens2)

MOCAPdata… data… MOCAP …myprecious… precious… …my



Time-Series and Matlab Time-series can be represented as vectors or arrays

Becoming Becomingsufficiently sufficiently familiar familiarwith withsomething something isisaa substitute substitute for for understanding understandingit. it.

– Fast vector manipulation • Most linear operations (eg euclidean distance, correlation) can be trivially vectorized

– Easy visualization – Many built-in functions – Specialized Toolboxes

•PART II: Time Series Matching Introduction



Basic Time-Series Matching Problem

Basic Data-Mining problem

Distance

Today’s databases are becoming too large. Search is difficult. How can we overcome this obstacle?

query D = 7.3

Basic structure of data-mining solution:

– Represent data in a new format

D = 10.2

Linear Scan:

– Search few data in the new representation

Objective: Compare the query with all sequences in DB and return the k most similar sequences to the query.

– Examine even fewer original data – Provide guarantees about the search results

D = 11.8

– Provide some type of data/result visualization D = 17

Database Databasewith withtime-series: time-series: – Medicalsequences sequences – Medical

– – Images, Images,etc etc

D = 22

Sequence SequenceLength:100-1000pts Length:100-1000pts DB DBSize: Size:11TByte TByte



Hierarchical Clustering

What other problems can we solve? Clustering: “Place time-series into ‘similar’ groups”



Very generic & powerful tool



Provides visual data grouping

Pairwise distances D1,1 D2,1

Classification: “To which group is a time-series most ‘similar’ to?” query ?

D M,N

? ?


1.

Merge objects with smallest distance

2.

Reevaluate distances

3.

Repeat process

Z = linkage(D); H = dendrogram(Z);


Partitional Clustering

K-Means Demo



Faster than hierarchical clustering



Typically provides suboptimal solutions (local minima)



Not good performance for high dimensions

1.4 1.2 1

K-Means Algorithm: 1. 2.

Initialize k clusters (k specified by user) randomly. Repeat until convergence 1. Assign each object to the nearest cluster center.

0.9

0.8

0.8

0.6

0.7

0.4

0.6

0.2

0.5

0

0.4

-0.2 0.3

-0.4

2. Re-estimate cluster centers.

0.2

-0.5 0

See: kmeans

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

0.5

1

1.5



K-Means Clustering for Time-Series

Classification Typically classification can be made easier if we have clustered the objects



So how is kMeans applied for Time-Series that are high-dimensional?



Perform kMeans on a compressed dimensionality

Class A 0.4

Q

Original sequences

Compressed sequences

0.2

0

Clustering space

-0.2

Project query in the new space and find its closest cluster

0.4

0.2

So, query Q is more similar to class B

-0.4

-0.6

- 0. 6

- 0. 4

- 0. 2

0

02 .

04 .

0 .6

08 .

Class B

0

-0.2

-0.4

-0.6

-0 .6

- 0. 4

- 0. 2

0

02 .

04 .

06 .

0 .8



Nearest Neighbor Classification

Example

We need not perform clustering before classification. We can classify an object based on the class majority of its nearest neighbors/matches.

Elfs Hobbits

10 9 8 7

h t g6 n e L5 r i a4 H

What do we need?

3

1. Define Similarity

2

2. Search fast

1 1

2

3

4

5

6

7

8

– Dimensionality Reduction (compress data)

9 10

Height



Notion of Similarity I 

Solution to any time-series problem, boils down to a proper definition of *similarity*

All Allmodels modelsare arewrong, wrong, but butsome someare areuseful… useful…

•PART II: Time Series Matching Similarity/Distance functions

Similarity is always subjective. (i.e. it depends on the application )



Notion of Similarity II

Metric and Non-metric Distance Functions

Similarity depends on the features we consider (i.e. how we will describe or compress the sequences)

Distance functions Metric

Non-Metric



Euclidean Distance



Time Warping



Correlation



LCSS

Properties Positivity: Positivity:d(x,y) d(x,y)≥≥00and andd(x,y)=0, d(x,y)=0,ififx=y x=y Symmetry: Symmetry:d(x,y) d(x,y)==d(y,x) d(y,x)

IfIfany anyof ofthese these is is not not obeyed obeyedthen thenthe thedistance distance is isaanon-metric non-metric

Triangle TriangleInequality: Inequality:d(x,z) d(x,z)≤≤d(x,y) d(x,y)++d(y,z) d(y,z)



Triangle Inequality

Triangle Inequality (Importance)

Triangle TriangleInequality: Inequality: d(x,z) d(x,z)≤≤d(x,y) d(x,y) ++d(y,z) d(y,z)

Triangle TriangleInequality: Inequality: d(x,z) d(x,z)≤≤d(x,y) d(x,y) ++ d(y,z) d(y,z) Assume:

Q

z Metric distance functions can exploit the triangle inequality to speed-up search

y

x

Intuitively, if: - x is similar to y and, - y is similar to z, then, - x is similar to z too.

d(Q,B) =150

A

d(Q,A) ≥ d(Q,B) – d(B,A)

B

So we don’t have to retrieve A from disk

d(Q,A) ≥ 150 – 20 = 130

C

A

B

C

A

0

20

110

B

20

0

90

C

110

90

0

80

100


Non-Metric Distance Functions

••Matching Matchingflexibility flexibility

Euclidean Distance 

Most widely used distance measure



Definition: L2 =

n

∑ (a[i] − b[i])

2

i =1

••Robustness Robustnessto tooutliers outliers ••Stretching Stretchingin intime/space time/space

Bat similar to batman

d(Q,bestMatch) = 20

Then, since d(A,B)=20


Man similar to bat??

and

••Support Supportfor fordifferent differentsizes/lengths sizes/lengths Batman similar to man

0

••Speeding-up Speeding-upsearch searchcan canbe be difficult difficult

20

40

60

uclidean distance distance L2 = sqrt(sum((a b).^2)); Euclidean b).^2)); % return E



Euclidean Distance (Vectorization)

Data Preprocessing (Baseline Removal)

Question: If I want to compare many sequences to each other do I have to use a for-loop? A

Answer: No, one can use the following equation to perform matrix computations only…

||A-B|| = sqrt ( ||A||2 + ||B||2 - 2*A.B )

B: DxN matrix Result is MxN matrix

B

M sequences

A: DxM matrix

A=

D h t g n e l f O

average value of A

average value of B

result D1,1 D2,1

…

D M,N

aa= .*b); ab=a'*b; aa=sum(a.*a); sum(a.*a); bb=sum(b bb= sum(b.*b); ab=a'*b; d = sqrt(repmat(aa',[1 size(bb,2)]) + repmat(bb,[size(aa,2) 1])

a = a – mean(a); mean(a);

- 2*ab ); 2*ab);



Data Preprocessing (Rescaling)

Dynamic Time-Warping (Motivation) Euclidean distance or warping cannot compensate for small distortions in time axis. A According to Euclidean distance B is more similar to A than to C

B C

Solution: Allow for compression & decompression in time

a = a ./ std(a); std(a);



Dynamic Time-Warping First used in speech recognition for recognizing words spoken at different speeds ---Maat--llaabb-------------------

Dynamic Time-Warping (how does it work?) Same idea can work equally well for generic time-series data

The intuition is that we copy an element multiple times so as to achieve a better matching

Euclidean Euclidean distance distance T1 T1==[1, [1,1, 1,2, 2,2] 2] dd==11

T2 T2==[1, [1,2, 2,2, 2,2] 2] One-to-one linear alignment ----Mat-lab--------------------------

Warping Warping distance distance T1 T1==[1, [1,1, 1,2, 2,2] 2] dd==00

T2 T2==[1, [1,2, 2,2, 2,2] 2] One-to-many non-linear alignment



Longest Common Subsequence (Implementation) Similar dynamic programming solution as DTW, but now we measure similarity not distance.

Distance Measure Comparison Dataset

Method

Time (sec)

Accuracy

Camera-Mouse

Euclidean

34

20%

DTW

237

80%

LCSS

210

100%

Euclidean

2.2

33%

DTW

9.1

44%

LCSS

8.2

46%

Euclidean

2.1

11%

DTW

9.3

15%

LCSS

8.3

31%

ASL

ASL+noise Can also be expressed as distance

LCSS offers enhanced robustness under noisy conditions



Distance Measure Comparison (Overview) Method Euclidean

Complexity

Elastic Matching

One-to-one Matching

Noise Robustness

O(n)

DTW

O(n* δ )

LCSS

O(n* δ )

•PART II: Time Series Matching Lower Bounding



Basic Time-Series Problem Revisited

Compression – Dimensionality Reduction Project all sequences into a new space, and search this space instead (egproject time- series from 100-D space to 2-D space)

Objective: Instead of comparing the query to the original sequences (Linear Scan/LS) , let’s compare the query to simplified versions of the DB time- series.

A 1 e r u t a e F

B C

Feature 2 query

query

This ThisDB DBcan cantypically typically fit fitin inmemory memory

One can also organize the low-dimensional points into a hierarchical ‘index’structure. In this tutorial we will not go over indexing techniques.

Question: When searching the original space it is guaranteed that we will find the best match. Does this hold (or under which circumstances) in the new compressed space?



Concept of Lower Bounding 

Generic Search using Lower Bounding

You can guarantee similar results to Linear Scan in the original dimensionality, as long as you provide a Lower Bounding (LB) function (in low dim) to the original distance (high dim.) GEMINI , GE neric M ultimedia IN dex In g

simplified DB

original DB

Answer Superset

Final Answer set

– So, for projection from high dim. (N) to low dim. (n): AÆa, BÆb etc Verify against original DB

D (a,b) <= D (A,B) DLB LB (a,b) <= Dtrue true(A,B) Projection onto X-axis

5

Α

C B

D

EF

4 1

C 3

F

B C

Β

1

Α

1

2

3

4

5

4

simplified query

Projection on some other axis

E

0 0

3

False alarm (not a problem)

D

2

1

2

2

D 3

EF 4

5

False dismissal (bad!)

5

query

“Find everything within range of 1 from A”



Lower Bounding Example

Lower Bounding Example query

sequences


Lower Bounding Example sequences

Lower Bounds

query

sequences


Lower Bounding Example sequences

Lower Bounds

True Distance

4.6399

4.6399

46.7790

37.9032

37.9032

108.8856

19.5174

19.5174

113.5873

72.1846

72.1846

104.5062

67.1436

67.1436

119.4087

78.0920

78.0920

120.0066

70.9273

70.9273

111.6011

63.7253

63.7253

119.0635

1.4121

1.4121

17.2540

BestSoFar



Lower Bounding the Euclidean distance

Fourier Decomposition “Everysignal signalcan can “Every berepresented representedas as be superpositionof of aasuperposition sinesand andcosines” cosines” sines (…alasnobody nobody (…alas believesme…) me…) believes

Decompose a time-series into sum of sine waves

There are many dimensionality reduction (compression ) techniques for time-series data. The following ones can be used to lower bound the Euclidean distance.

DFT: IDFT:

0

20 40 60 80 100 120

DFT

0

20 40 6 0 80 1 00 120

0

20 4 0 6 0 8 0 1 00 12 0

DWT

0

SVD

2 0 4 0 6 0 8 0 1 00 12 0

0

2 0 4 0 6 0 8 0 1 00 12 0

APCA

PAA

0

2 0 4 0 6 0 8 0 1 00 12 0

PLA

Figure by Eamonn Keogh, ‘Time-Series Tutorial’



Fourier Decomposition Decompose a time-series into sum of sine waves DFT:

X(f)

x(n)

-0.3633

-0.4446

-0.6280 + 0.2709i

-0.9864

-0.4929 + 0.0399i

-0.3254

-1.0143 + 0.9520i

-0.6938

0.7200 -1.0571i

-0.1086

-0.0411 + 0.1674i

-0.3470

-0.5120 -0.3572i

IDFT:

decomposition fa = fft(a); fft(a); % Fourier decomposition fa(5:end) = 0; % keep first 5 coefficients (low frequencies) reconstr = real(ifft(fa )); % reconstruct signal

Fourier Decomposition How much space we gain by compressing random walk data?

0.5849

0.9860 + 0.8043i

1.5927

-0.3680 -0.1296i

-0.9430

-0.0517 -0.0830i

-0.3037

-0.9158 + 0.4481i

-0.7805

1.1212 -0.6795i

-0.1953

0.2667 + 0.1100i

-0.3037

0.2667 -0.1100i

0.2381

1.1212 + 0.6795i

2.8389

-0.9158 -0.4481i

-0.7046

-0.0517 + 0.0830i

-0.5529

-0.3680 + 0.1296i

-0.6721

0.9860 -0.8043i

0.1189

-0.5120 + 0.3572i

0.2706

-0.0411 -0.1674i

-0.0003

0.7200 + 1.0571i

1.3976

-1.0143 -0.9520i

-0.4987

-0.4929 -0.0399i

-0.2387

-0.6280 -0.2709i

-0.7588

Reconstruction using 1coefficients

5

0

-5 50



1 coeff > 60% of energy




100

150

200

250

Life is complex, it has both real and imaginary parts.



Fourier Decomposition


How much space we gain by compressing random walk data?




5

5

0

0

-5

-5 50

100

150

200

250

50













100

150

200

250






How much space we gain by compressing random walk data? EnergyPercentage

Error 1

1500

0.95 0.9

Reconstruction using 20coefficients 1000

0.85

5

0.8 0.75

500

0

0.7 -5

0.65 50

100

150

200

0

250

20

40

60 80 Coefficients

100

120














20

40

60 80 Coefficients

100

120




Which coefficients are important?

Which coefficients are important?

– We can measure the ‘energy’ of each coefficient

– We can measure the ‘energy’ of each coefficient

– Energy = Real(X(f k ))2 + Imag(X(f k ))2

– Energy = Real(X(f k ))2 + Imag(X(f k ))2 Most of data-mining research uses first k coefficients: 

Good for random walk signals (eg stock market)



Easy to ‘index’



Not good for general signals

Usage of the coefficients with highest energy: 

Good for all types of signals



Believed to be difficult to index



CAN be indexed using metric trees

fa = fft(a); decomposition fft(a); % Fourier decomposition N = length(a ); % how many? many? fa = fa(1:ceil(N/2)); % keep first half only mag = 2*abs(fa).^2; % calculate energy energy



Code for Reconstructed Sequence a = load('randomWalk.dat'); a = (a(a mean(a))/std(a); mean(a))/std(a);

X(f) 0 -0.6280 + 0.2709i

% zz-normalization

keep

-0.4929 + 0.0399i -1.0143 + 0.9520i 0.7200 -1.0571i

fa = fft(a);

-0.0411 + 0.1674i

maxInd = ceil(length(a)/2); N = length(a);

% until the middle

-0.5120 -0.3572i 0.9860 + 0.8043i

Code for Plotting the Error a = load('randomWalk.dat'); a = (a(a mean(a))/std(a); mean(a))/std(a); fa = fft(a); maxInd = ceil(length(a)/2); N = length(a); energy = zeros(maxInd -1, 1); E = sum(a.^2);

% zz-normalization This is the same

% until the middle % energy of a

-0.3680 -0.1296i

energy = zeros(maxIndzeros(maxInd-1, 1); E = sum(a.^2);

-0.0517 -0.0830i

% energy of a

1.1212 -0.6795i

for ind=2:maxInd, fa_N = fa; fa_N(ind+1:N -ind+1) = 0; r = real(ifft(fa_N));

end

-0.9158 + 0.4481i

Ignore % copy fourier % zero out unused % reconstruction

plot(r, 'r','LineWidth',2); hold on; plot(a,'k'); title(['Reconstruction using ' num2str(indnum2str(ind-1) 'coefficients']); set(gca,'plotboxaspectratio', [3 1 1]); axis tight pause; % wait for key cla; % clear axis keep

0.2667 + 0.1100i 0.2667 -0.1100i 1.1212 + 0.6795i -0.9158 -0.4481i -0.0517 + 0.0830i -0.3680 + 0.1296i 0.9860 -0.8043i -0.5120 + 0.3572i

for ind=2:maxInd, fa_N = fa; fa_N(ind+1:N -ind+1) = 0; r = real(ifft(fa_N));

% % %

copy fourier zero out unused reconstruction

energy(ind -1) = sum(r.^2); % energy of reconstruction error(ind -1) = sum(abs(r -a).^2); % error end E = ones(maxInd -1, 1)*E; error = E - energy; ratio = energy ./ E;

-0.0411 -0.1674i 0.7200 + 1.0571i -1.0143 -0.9520i -0.4929 -0.0399i -0.6280 -0.2709i

subplot(1,2,1); % left plot plot([1:maxInd -1], error, 'r', 'LineWidth',1.5); subplot(1,2,2); % right plot plot([1:maxInd -1], ratio, 'b', 'LineWidth',1.5);



Lower Bounding using Fourier coefficients

Lower Bounding using Fourier coefficients -Example

Parseval’s Theorem states that energy in the frequency domain equals the energy in the time domain:

x y Euclidean distance

or, that

If we just keep some of the coefficients, their sum of squares always underestimates (ie lower bounds) the Euclidean distance:

Note the normalization x = cumsum(randn(100,1)); y = cumsum(randn(100,1)); euclid_Time = sqrt(sum((xsqrt(sum((x-y).^2));

120.9051

fx = fft(x)/sqrt(length(x)); fft(x)/sqrt(length(x)); fy = fft(y)/sqrt(length(x)); fft(y)/sqrt(length(x)); euclid_Freq = sqrt(sum(abs(fx - fy).^2));


O(nlogn) O(nlogn)complexity complexity



Tried Triedand andtested tested



Hardware Hardwareimplementations implementations



Many Manyapplications: applications:

– – compression compression – smoothing – smoothing

Wavelets – Why exist?



Not Notgood goodapproximation approximationfor for bursty signals bursty signals



Not Notgood goodapproximation approximationfor for signals signalswith withflat flatand andbusy busy sections sections (requires (requiresmany manycoefficients) coefficients)



Similar concept with Fourier decomposition



Fourier coefficients represent global contributions, wavelets are localized

Fourier is good for smooth, random walk data, but not for bursty data or flat data

– – periodicity periodicitydetection detection



Wavelets (Haar) - Intuition

Wavelets in Matlab



Wavelet coefficients, still represent an inner product (projection) of the signal with some basis functions.



These functions have lengths that are powers of two (full sequence length, half, quarter etc)

Specialized Matlab interface for wavelets

An arithmetic example

c-d00

X = [9,7,3,5] c+d00 D

etc

Haar coefficients: {c, d00, d10, d11,…}

Haar = [6,2,1,-1]

c = 6 = (9+7+3+5)/4 c + d00 = 6+2 = 8 = (9+7)/2 c - d00 = 6-2 = 4 = (3+5)/2 etc

See also :wavemenu

120.9051





Keeping 10 coefficients the distance is: 115.5556 < 120.9051

See also:wavemenu



Code for Haar Wavelets

PAA (Piecewise Aggregate Approximation)

a = load('randomWalk.dat'); a = (a% z(a mean(a))/std(a); mean(a))/std(a); z-normalization maxlevels = wmaxlev(length(a),'haar'); [Ca, La] = wavedec(a,maxlevels,'haar');



Represent time-series as a sequence of segments



Essentially a projection of the Haar coefficients in time

also featured as Piecewise Constant Approximation

% Plot coefficients and MRA for level = 1:maxlevels cla; subplot(2,1,1); plot(detcoef(Ca,La,level)); axis tight; title(sprintf('Wavelet coefficients – Level %d',level)); subplot(2,1,2); plot(wrcoef('d',Ca,La,'haar',level)); plot(wrcoef('d',Ca,La,'haar',level)); axis tight; title(sprintf('MRA – Level %d',level)); pause; end

Reconstruction using 1coefficients 2 1 0 -1 -2 50

100

150

200

250

% TopTop-20 coefficient reconstruction [Ca_sorted, Ca_sortind] = sort(Ca); Ca_top20 = Ca; Ca_top20(Ca_sortind(1:end Ca_top20(Ca_sortind(1:end -19)) = 0; a_top20 = waverec(Ca_top20,La,'haar'); figure; hold on; plot(a, 'b'); plot(a_top20, 'r');



PAA (Piecewise Aggregate Approximation) also featured as Piecewise Constant Approximation
















2

2

1

1

0

0

-1

-1

-2

-2 50

100

150

200

250

50


100

150

200

250


















2

2

1

1

0

0

-1

-1

-2

-2 50

100

150

200

250

50

100

150

200

250



PAA (Piecewise Aggregate Approximation)

PAA Matlab Code

also featured as Piecewise Constant Approximation







function data = paa(s, numCoeff) % PAA(s, numcoeff) % s: sequence vector (Nx1 or Nx1) % numCoeff: number of PAA segments % data: PAA sequence (Nx1) N = length(s); segLen = N/numCoeff;

Reconstruction using 32coefficients 2

% length of sequence % assume it's integer

sN = reshape(s, segLen, numCoeff); avg = mean(sN); data = repmat(avg, segLen, 1); data = data(:);

1 0

% % % %

break in segments average segments expand segments make column

-1 -2 50

100

150

200

s

250

1


2

3

4

5

6

7

PAA Matlab Code

function data = paa(s, numCoeff) % PAA(s, numcoeff) % s: sequence vector (Nx1 or Nx1) % numCoeff: number of PAA segments % data: PAA sequence (Nx1)




s

1

2

3

4

5

6

7

% % % %


numCoeff

8

N=8 segLen = 2

N = length(s); segLen = N/numCoeff;

% length of sequence % assume it's integer 2

4

s

1

2

3

4

1

3

5

7

2

4

6

8


5

6

7

PAA Matlab Code


function data = paa(s, numCoeff) % PAA(s, numcoeff) % s: sequence vector (1xN) % numCoeff: number of PAA segments % data: PAA sequence (1xN)



sN avg

1

2

3

4

1

3

5

7

2

4

6

8

1.5

3.5

5.5

7.5

5

6

7

% % % %


numCoeff

8

4


PAA Matlab Code

s

N=8 segLen= 2

4


sN


4


PAA Matlab Code


numCoeff

8

% % % %

8

N=8 segLen = 2


numCoeff



sN = reshape(s, segLen, numCoeff); avg = mean(sN); 2 data = repmat(avg, segLen, 1); data = data(:) ’;

4

s sN avg

1

2

3

4

1

3

5

7

2

4

6

8

1.5

3.5

5.5

7.5

5

6

7

% % % %

break in segments average segments expand segments make row

numCoeff

8

data

N=8 segLen= 2

1.5

3.5

5.5

7.5

1.5

3.5

5.5

7.5

4



PAA Matlab Code

APCA (Adaptive Piecewise Constant Approximation) PAA

function data = paa(s, numCoeff) % PAA(s, numcoeff) % s: sequence vector (1xN) % numCoeff: number of PAA segments % data: PAA sequence (1xN) N = length(s); segLen = N/numCoeff;

Segments of equal size


sN = reshape(s, segLen, numCoeff); avg = mean(sN); data = repmat(avg, segLen, 1); data = data(:) ’;

% % % %

N=8 segLen = 2

break in segments average segments expand segments make row APCA

s sN avg

1

2

3

4

1

3

5

7

2

4

6

8

1.5

3.5

5.5

7.5

5

6

7

numCoeff

8

data data

1.5

3.5

5.5

7.5

1.5

3.5

5.5

7.5

1.5

1.5

3.5

3.5

Segments of variable size

4

5.5

5.5

7.5



Not all haar/PAA coefficients are equally important



Intuition: Keep ones with the highest energy



Segments of variable length



APCA is good for bursty signals



PAA requires 1 number per segment, APCA requires 2: [value, length] E.g. 10 bits for a sequence of 1024 points

7.5



Wavelet Decomposition

Piecewise Linear Approximation (PLA) 

Approximate a sequence with multiple linear segments



First such algorithms appeared in cartography for map approximation



Many implementations



O(n) O(n)complexity complexity



Hierarchical Hierarchicalstructure structure



Progressive Progressivetransmission transmission



Better Betterlocalization localization

– Greedy Bottom-Up



Good Goodfor forbursty burstysignals signals

– Greedy Top-down



Many Manyapplications: applications:

– Genetic, etc



Most Mostdata-mining data-miningresearch research still stillutilizes utilizes Haar Haarwavelets wavelets because becauseof oftheir theirsimplicity. simplicity.

– – compression compression – – periodicity periodicitydetection detection

– Optimal



You can find a bottom-up implementation here:

– http://www.cs.ucr.edu/~eamonn/TSDMA/time_series_toolbox/



Piecewise Linear Approximation (PLA)










































O(nlogn) O(nlogn)complexity complexityfor for “bottom “bottomup” up”algorithm algorithm



Incremental Incrementalcomputation computation possible possible



Provable Provableerror errorbounds bounds



Applications Applicationsfor: for:



Visually Visuallynot notvery verysmooth smoothor or pleasing. pleasing.

– – Image Image/ /signal signal simplification simplification – – Trend Trenddetection detection



Singular Value Decomposition (SVD)

Singular Value Decomposition (SVD)



SVD attempts to find the ‘optimal’ basis for describing a set of multidimensional points



Objective: Find the axis (‘directions’) that describe better the data variance



Each time-series is essentially a multidimensional point



Objective: Find the ‘eigenwaves’ (basis) whose linear combination describes best the sequences. Eigenwaves are data-dependent.

AMxn = UMxr *Σ rxr * VTnxr

eigenwave 0

Factoring of data array into 3 matrices x

eigenwave 1

x

each of length n n eigenwave 3

eigenwave 4

y We need 2 numbers (x,y) for every point

y Now we can describe each point with 1 number, their projection on the line

New axis and position of points (after projection and rotation)

A linear combination of the eigenwavescan produce any sequence in the database

s e c n e u q e s M

[U,S,V] = svd(A) svd(A)

…



Code for SVD / PCA

Singular Value Decomposition

A = cumsum(randn(100,10)); % zz-normalization A = (A -repmat(mean(A),size(A,1),1))./repmat(std(A),size(A,1),1); [U,S,V] = svd(A,0); % Plot relative energy figure; plot(cumsum(diag(S).^2)/norm(diag(S))^2); set(gca, 'YLim', [0 1]); pause; % TopTop-3 eigenvector reconstruction A_top3 = U(:,1:3)*S(1:3,1:3)*V(:,1:3)'; % Plot original and reconstruction figure; for i = 1:10 cla; subplot(2,1,1); plot(A(:,i)); title('Original'); axis tight; subplot(2,1,2); plot(A_top3(:,i)); title('Reconstruction'); axis tight; pause; end



Optimal Optimaldimensionality dimensionality reduction reductionin inEuclidean Euclidean distance distancesense sense



SVD SVDis isaavery verypowerful powerfultool tool in inmany manydomains: domains:



Cannot Cannotbe beapplied appliedfor forjust just one onesequence. sequence.AAset setof of sequences is required. sequences is required.



Addition Additionof ofaasequence sequencein in database databaserequires requires recomputation recomputation



Very Verycostly costlyto tocompute. compute. 2 2 Time: Time:min{ min{O(M O(M 2n), n),O(Mn O(Mn2)})} Space: Space:O(Mn) O(Mn)

– – Websearch Websearch(PageRank) (PageRank)

MMsequences sequencesof oflength lengthnn



Symbolic Approximation

Symbolic Approximations



Assign a different symbol based on range of values



Find ranges either from data histogram or uniformly

c c

c b

b a

0

20

b

a 40

60

80

100

Linear Linearcomplexity complexity



After After‘symbolization’ ‘symbolization’many many tools toolsfrom frombioinformatics bioinformatics can canbe beused used



Number Numberof ofregions regions (alphabet (alphabetlength) length)can canaffect affect the thequality qualityof ofresult result

– – Markov Markovmodels models – – Suffix-Trees, Suffix-Trees,etc etc

120

baabccbc 



You can find an implementation here:

– http://www.ise.gmu.edu/~jessica/sax.htm


Multidimensional Time-Series 

Catching momentum lately



Applications for mobile trajectories, sensor networks, epidemiology , etc



Let’s see how to approximate 2D trajectories with Minimum Bounding Rectangles


Ari,are areyou yousure surethe the Ari, worldisisnot not1D? 1D? world

Aristotle

Multidimensional MBRs Find Bounding rectangles that completely contain a trajectory given some optimization criteria (eg minimize volume)

On my income tax 1040 it says "Check this box if you are blind." I wanted to put a check mark about three inches away. - Tom Lehrer



Comparison of different Dim. Reduction Techniques

So which dimensionality reduction is the best? Fourierisis Fourier good… good…

PAA! PAA!

1993

APCAisis APCA better better thanPAA! PAA! than

2001

2000

Chebyshev Chebyshev better isisbetter thanAPCA! APCA! than

2004

The The futureisis future symbolic! symbolic!

2005

Absence of proof is no proof of absence. - Michael Crichton



Comparisons Lets see how tight the lower bounds are for a variety on 65 datasets Average Lower Bound

A. No approach is better on all datasets B. Best coeff. techniques can offer tighter bounds

Median Lower Bound

•PART II: Time Series Matching Lower Bounding the DTW and LCSS

C. Choice of compression depends on application Note: similar results also reported by Keogh in SIGKDD02



Lower Bounding the Dynamic Time Warping


Recent approaches use the Minimum Bounding Envelope for bounding the DTW

LB by Keogh approximate MBE and sequence using MBRs

– Create Minimum Bounding Envelope (MBE) of query Q – Calculate distance between MBE of Q and any sequence A

LB = 13.84

)δ ,A) < DTW(Q,A) – One can show that: D(MBE(Q D(MBE(Q Q

A

LB = sqrt(sum([[A > U].* [A-U]; [A < L].* [L-A]].^2)); One Matlab command! U

δ

MBE(Q)

LB by Zhu and Shasha approximate MBE and sequence using PAA

A

Q

L

However, this representation is uncompressed. Both MBE and the DB sequence can be compressed using any of the previously mentioned techniques.

LB = 25.41 Q A




Time Comparisons

An even tighter lower bound can be achieved by ‘warping’ the MBE approximation against any other compressed signal.

We will use DTW (and the corresponding LBs) for recognition of hand-written digits/shapes.

LB_Warp = 29.05

Lower Bounding approaches for DTW, will typically yield at least an order of magnitude speed improvement compared to the naïve approach.

Accuracy: Using DTW we can achieve recognition above 90%. Running Time: runTime LB_Warp < runTimeLB_Zhu < runTimeLB-Keogh Pruning Power: For some queries LB_Warp can examine up to 65 time fewer sequences

Let’s compare the 3 LB approaches:



Upper Bounding the LCSS

LCSS Application – Image Handwriting

Since LCSS measures similarity and similarity is the inverse of distance, to speed up LCSS we need to upper bound it.



Library of Congress has 54 million manuscripts (20TB of text)



Increasing interest for automatic transcribing

Word annotation: 1.1.Extract Extractwords wordsfrom fromdocument document 2.2.Extract Extractimage imagefeatures features 3.3.Annotate Annotateaasubset subsetof ofwords words 4.4.Classify Classifyremaining remainingwords words

LCSS(MBE >= LCSS(Q,A) Q,A) LCSS(MBE Q,A) >= LCSS(Q,A) Query

Indexed Sequence

1

0.8

Sim.=50/77 = 0.64

e u l a0.6 V e r u t a0.4 e F

0.2

0

50

100

150

200

250

300

350

Column

Features:

44 points

+

6 points George Washington Manuscript


- Black pixels / column - Ink-paper transitions/ col , etc


LCSS Application – Image Handwriting Utilized 2D time-series (2 features) Returned 3-Nearest Neighbors of following words Classification accuracy > 70%

•PART II: Time Series Analysis Test Case and Structural Similarity Measures

400



Analyzing Time-Series Weblogs

Weblog Data Representation Record aggregate information, eg, number of requests per day for each We can keyword Query: Spiderman

“PKDD 2005”

May 2002. Spiderman 1 was released in theaters

s t s e u q e R

“Porto”

Jan

Privacy preserving

Nov Dec

Google Zeitgeist

Finding similar patterns in query logs

We can find useful patterns and correlation in the user demand patterns which can be useful for:

We can find useful patterns and correlation in the user demand patterns which can be useful for:



Search engine optimization





Recommendations



Advertisement pricing (e.g. keyword more expensive at the popular months)

Query: xbox

s t s e u q e R

Query: ps2 Apr May Jun

Aug Sep Okt


Finding similar patterns in query logs

Feb Mar

Jul




Jan

Apr May Jun

Capture trends and periodicities

Weblog of user requests over time

“Priceline”



Feb Mar



Jul

Aug Sep Okt

Nov Dec



Search engine optimization Recommendations Advertisement pricing (e.g. keyword more expensive at the popular months)

s t s e u q e R

Query: elvis Jan

Feb Mar

Apr May Jun

Jul

Nov Dec

th Burst on Aug. 16 Death Anniversary of Elvis

Game consoles are more popular closer to Christmas


Aug Sep Okt


Matching of Weblog data

First Fourier Coefficients vs Best Fourier Coefficients

Use Euclidean distance to match time-series. But which dimensionality reduction technique to use? Let’s look at the data:

Query “Bach” 1 year span

The data is smooth and highly periodic, so we can use Fourier decomposition. Instead of using the first Fourier coefficients we can use the best ones instead. Let’s see how the approximation will look:

Query “stock market”

Using the best coefficients, provides a very high quality approximation of the original time-series



Matching results I

Matching results II

Query = “Lance Armstrong”

2000

2001

Query = “Christmas”

2002

2000

2001

2002

Knn4: Christmas coloring books

LeTour 0

Knn8: Christmas baking 2000

2001

2002

Knn12: Christmas clipart Tour De France Knn20: Santa Letters 0

2000

2001

2002



Finding Structural Matches

Periodic Matching

The Euclidean distance cannot distill all the potentially useful information in the weblog data. 

Some data are periodic , while other are bursty . We will attempt to provide similarity measures that are based on periodicity and burstiness.

Frequency

F ( x), F ( y)

Ignore Phase/ Keep important components

Calculate Distance

arg max || F ( x) ||, F ( x + ) k arg max || F ( y ) ||, F ( y + ) k

cinema

D1 =|| F ( x + ) − F ( y + ) ||

Periodogram

D2 =|| F ( x + ) ⋅ F ( y + ) ||

Query “cinema”. Weakly periodicity. Peak of period every Friday.

stock

easter

Query “Elvis”. Burst in demand on 16 th August. Death anniversary of Elvis Presley


Matching Results with Periodic Measure Now we can discover more flexible matches. We observe a clear separation between seasonal and periodic sequences.

0

5

10

15

20

25

30

35

40

45

50

0

5

10

15

20

25

30

35

40

45

50

christmas


Matching Results with Periodic Measure Compute pairwise periodic distances and do a mapping of the sequences on 2D using Multi-dimensional scaling (MDS).

Time Series Analysis With Matlab Tutorials

Recommend Documents