Processing Independent Component Analysis Class 8. 23 Sep 2013 - - PowerPoint PPT Presentation

Machine Learning for Signal Processing Independent Component Analysis

Class 8. 23 Sep 2013 Instructor: Bhiksha Raj

23 Sep 2013 11755/18797 1

Correlation vs. Causation

• The consumption of burgers has gone up

steadily in the past decade

• In the same period, the penguin population of

Antarctica has gone down

23 Sep 2013 11755/18797 2

The concept of correlation

• Two variables are correlated if knowing the

value of one gives you information about the expected value of the other

23 Sep 2013 11755/18797 3

Burger consumption Penguin population Time

The statistical concept of correlatedness

• Two variables X and Y are correlated if If

knowing X gives you an expected value of Y

• X and Y are uncorrelated if knowing X tells you

nothing about the expected value of Y

– Although it could give you other information – How?

23 Sep 2013 11755/18797 4

A brief review of basic probability

• Uncorrelated: Two random variables X and Y are

uncorrelated iff:

– The average value of the product of the variables equals the product of their individual averages

• Setup: Each draw produces one instance of X and one

instance of Y

– I.e one instance of (X,Y)

• E[XY] = E[X]E[Y]
• The average value of X is the same regardless of the value
• f Y

23 Sep 2013 11755/18797 5

Uncorrelatedness

• Which of the above represent uncorrelated RVs?

23 Sep 2013 11755/18797 6

The statistical concept of Independence

• Two variables X and Y are dependent if If

knowing X gives you any information about Y

• X and Y are independent if knowing X tells you

nothing at all of Y

23 Sep 2013 11755/18797 7

A brief review of basic probability

• Independence: Two random variables X and Y are

independent iff:

– Their joint probability equals the product of their individual probabilities

• P(X,Y) = P(X)P(Y)
• Independence implies uncorrelatedness

– The average value of X is the same regardless of the value of Y

• E[X|Y] = E[X]

– But not the other way

23 Sep 2013 11755/18797 8

A brief review of basic probability

• Independence: Two random variables X and Y

are independent iff:

• The average value of any function of X is the

same regardless of the value of Y

– Or any function of Y

• E[f(X)g(Y)] = E[f(X)] E[g(Y)] for all f(), g()

23 Sep 2013 11755/18797 9

Independence

• Which of the above represent independent RVs?
• Which represent uncorrelated RVs?

23 Sep 2013 11755/18797 10

A brief review of basic probability

• The expected value of an odd function of an

RV is 0 if

– The RV is 0 mean – The PDF is of the RV is symmetric around 0

• E[f(X)] = 0 if f(X) is odd symmetric

23 Sep 2013 11755/18797 11

y = f(x) p(x)

A brief review of basic info. theory

• Entropy: The minimum average number of bits

to transmit to convey a symbol

• Joint entropy: The minimum average number of

bits to convey sets (pairs here) of symbols

23 Sep 2013 11755/18797 12

T(all), M(ed), S(hort)… T, M, S… M F F M..



 

X

X P X P X H )] ( log )[ ( ) (



 

Y X

Y X P Y X P Y X H

,

)] , ( log )[ , ( ) , ( X Y

A brief review of basic info. theory

• Conditional Entropy: The minimum average

number of bits to transmit to convey a symbol X, after symbol Y has already been conveyed

– Averaged over all values of X and Y

23 Sep 2013 11755/18797 13

T, M, S… M F F M.. X Y

  

   

Y X Y X

Y X P Y X P Y X P Y X P Y P Y X H

,

)] | ( log )[ , ( )] | ( log )[ | ( ) ( ) | (

A brief review of basic info. theory

• Conditional entropy of X = H(X) if X is

independent of Y

• Joint entropy of X and Y is the sum of the

entropies of X and Y if they are independent

23 Sep 2013 11755/18797 14

) ( )] ( log )[ ( ) ( )] | ( log )[ | ( ) ( ) | ( X H X P X P Y P Y X P Y X P Y P Y X H

Y X Y X

    

     

   

Y X Y X

Y P X P Y X P Y X P Y X P Y X H

, ,

)] ( ) ( log )[ , ( )] , ( log )[ , ( ) , ( ) ( ) ( ) ( log ) , ( ) ( log ) , (

, ,

Y H X H Y P Y X P X P Y X P

Y X Y X

     



Onward..

23 Sep 2013 11755/18797 15

Projection: multiple notes

23 Sep 2013 11755/18797 16

 P = W (WTW)-1 WT  Projected Spectrogram = PM

M = W =

We’re actually computing a score

23 Sep 2013 11755/18797 17

 M ~ WH  H = pinv(W)M

M = W = H = ?

So what are we doing here?

• M ~ WH is an approximation
• Given W, estimate H to minimize error
• Must ideally find transcription of given notes

 



   

i j ij ij F 2 2

) ( min arg || || min arg H W M H W M H

H H

W = H = ? ?

23 Sep 2013 11755/18797 18

How about the other way?

23 Sep 2013 11755/18797 19

 M ~ WH W = Mpinv(H) U = WH

M = W =

? ?

H = U =

Going the other way..

• M ~ WH is an approximation
• Given H, estimate W to minimize error
• Must ideally find the notes corresponding to the

transcription

 



   

i j ij ij F 2 2

) ( min arg || || min arg H W M H W M W

H W

W =? H ?

23 Sep 2013 11755/18797 20

When both parameters are unknown

• Must estimate both H and W to best

approximate M

• Ideally, must learn both the notes and their

transcription!

W =? H = ? approx(M) = ?

23 Sep 2013 11755/18797 21

A least squares solution

• Unconstrained

– For any W, H that minimizes the error, W’=WA, H’=A-1H also minimizes the error for any invertible A

• Too many solutions

2 ,

|| || min arg ,

F

H W M H W

H W

 

H

23 Sep 2013 11755/18797 22

A constrained least squares solution

• For our problem, lets consider the “truth”..
• When one note occurs, the other does not

– hi

Thj = 0 for all i != j

• The rows of H are uncorrelated

2 ,

|| || min arg ,

F

H W M H W

H W

 

H

A least squares solution

• Assume: HHT = I

– Normalizing all rows of H to length 1

• pinv(H) = HT
• Projecting M onto H

– W = M pinv(H) = MHT – WH = M HTH

H

2 ,

|| || min arg ,

F

H W M H W

H W

 

2

|| || min arg

F TH

H M M H

H

 

Constraint: Rank(H) = 4

Finding the notes

• Add the constraint: HHT = I
• The solution is obtained through Eigen

decomposition

• Note: we are considering the correlation of

MT

23 Oct 2012 11755/18797 26

  H H M ) (

T

n Correlatio

 

T F T

H H H H M M H

H

   

2

|| || min arg

So how does that work?

• There are 12 notes in the segment, hence we

try to estimate 12 notes..

23 Oct 2012 11755/18797 28

So how does that work?

• The scores of the first three “notes” and their

contributions

23 Oct 2012 11755/18797 29

Finding the notes

• Can find W instead of H
• Assume the columns of W are orthogonal
• This results in the more conventional Eigen

decomposition

23 Oct 2012 11755/18797 30

  W W M) ( n Correlatio

2

|| || min arg

F T

M W W M W

W

 

So how does that work?

• There are 12 notes in the segment, hence we

try to estimate 12 notes..

• Results are not good again

23 Oct 2012 11755/18797 31

Our notes are not orthogonal

• Overlapping frequencies
• Note occur concurrently

– Harmonica continues to resonate to previous note

• More generally, simple orthogonality will not give

us the desired solution

23 Sep 2013 11755/18797 32

Eigendecomposition and SVD

• Matrix M can be decomposed as M = USVT
• When we assume the scores are orthogonal, we get

H = VT, W = US

• When we assume the notes are orthogonal, we get

W = U, H = SVT

• In either case the results are the same

– The notes are orthogonal and so are the scores – Not good in our problem

23 Sep 2013 11755/18797 34

T

USV M  WH M 

Orthogonality

• In any least-squared error decomposition

M=WH, if the columns of W are orthogonal, the rows of H will also be orthogonal

• Sometimes mere orthogonality is not enough

23 Sep 2013 11755/18797 35

WH M 

What else can we look for?

• Assume: The “transcription” of one note does

not depend on what else is playing

– Or, in a multi-instrument piece, instruments are playing independently of one another

• Not strictly true, but still..

23 Sep 2013 11755/18797 36

Formulating it with Independence

• Impose statistical independence constraints
• n decomposition

23 Sep 2013 11755/18797 37

) . . . . ( || || min arg ,

2 ,

t independen are H

• f

rows

F

    H W M H W

H W

Changing problems for a bit

• Two people speak simultaneously
• Recorded by two microphones
• Each recorded signal is a mixture of both signals

23 Sep 2013 11755/18797 38

) ( ) ( ) (

2 12 1 11 1

t h w t h w t m   ) ( ) ( ) (

2 22 1 21 2

t h w t h w t m  

) (

1 t

h ) (

2 t

h

A Separation Problem

• M = WH

– M = “mixed” signal – W = “notes” – H = “transcription”

• Separation challenge: Given only M estimate H
• Identical to the problem of “finding notes”

23 Sep 2013 11755/18797 39

=

M H W w11 w12 w21 w22 Signal at mic 1 Signal at mic 2 Signal from speaker 1 Signal from speaker 2

Imposing Statistical Constraints

• M = WH
• Given only M estimate H
• H = W-1M = AM
• Only known constraint: The rows of H are

independent

• Estimate A such that the components of AM are

statistically independent

– A is the unmixing matrix

23 Sep 2013 11755/18797 40

=

M H W w11 w12 w21 w22

Statistical Independence

• M = WH H = AM
• Emulating independence

– Compute W (or A) and H such that H has statistical characteristics that are observed in statistically independent variables

• Enforcing independence

– Compute W and H such that the components of M are independent

23 Sep 2013 11755/18797 41

Remember this form

Emulating Independence

• The rows of H are uncorrelated

– E[hihj] = E[hi]E[hj] – hi and hj are the ith and jth components of any vector in H

• The fourth order moments are independent

– E[hihjhkhl] = E[hi]E[hj]E[hk]E[hl] – E[hi

2hjhk] = E[hi 2]E[hj]E[hk]

– E[hi

2hj 2] = E[hi 2]E[hj 2]

– Etc.

23 Sep 2013 11755/18797 42

H

Zero Mean

• Usual to assume zero mean processes

– Otherwise, some of the math doesn’t work well

• M = WH H = AM
• If mean(M) = 0 => mean(H) = 0

– E[H] = A.E[M] = A0 = 0 – First step of ICA: Set the mean of M to 0 – mi are the columns of M

23 Sep 2013 11755/18797 43





i i

cols m M

m

) ( 1 

i

i i

  

m

m m 

Emulating Independence..

• Independence  Uncorrelatedness
• Estimate a C such that CM is uncorrelated
• X = CM

– E[xixj] = E[xi]E[xj] = dij [since M is now “centered”] – XXT = I

• In reality, we only want this to be a diagonal matrix, but we’ll

make it identity

23 Sep 2013 11755/18797 44

H H’

=

Diagonal + rank1 matrix H=AM H=BCM A=BC

Decorrelating

• X = CM
• XXT = I
• Eigen decomposition MMT= ESET
• Let C = S-1/2ET

– X = S-1/2ETM – XXT = WMMTWT = S-1/2ET ESETES-1/2 = I

23 Sep 2013 11755/18797 45

H H’

=

Diagonal + rank1 matrix H=AM H=BCM A=BC

Decorrelating

• X = CM
• XXT = I
• Eigen decomposition MMT= ESET
• Let C = S-1/2ET

– X = S-1/2ETM – WMMTWT = S-1/2ET ESETES-1/2 = I

• X is called the whitened version of M

– The process of decorrelating M is called whitening – C is the whitening matrix

23 Sep 2013 11755/18797 46

H H’

=

Diagonal + rank1 matrix H=AM H=BCM A=BC

Uncorrelated != Independent

• Whitening merely ensures that the resulting shat

signals are uncorrelated, i.e. E[xixj] = 0 if i != j

• This does not ensure higher order moments are also

decoupled, e.g. it does not ensure that E[xi

2xj 2] = E[xi 2]E [xj 2]

• This is one of the signatures of independent RVs
• Lets explicitly decouple the fourth order moments

47 23 Sep 2013 11755/18797

Decorrelating

• X = CM
• XXT = I
• Will multiplying X by B re-correlate the components?
• Not if B is unitary

– BBT = BTB = I

• HHT = BXXTBT = BBT = I
• So we want to find a unitary matrix

– Since the rows of H are uncorrelated

• Because they are independent

23 Sep 2013 11755/18797 48

H H’

=

Diagonal + rank1 matrix H=AM H=BX A=BC

ICA: Freeing Fourth Moments

• We have E[xi xj] = 0 if i != j

– Already been decorrelated

• A=BC, H = BCM, X = CM,  H = BX
• The fourth moments of H have the form:

E[hi hj hk hl]

• If the rows of H were independent

E[hi hj hk hl] = E[hi] E[hj] E[hk] E[hl]

• Solution: Compute B such that the fourth moments of H = BX

are decoupled

– While ensuring that B is Unitary

49 23 Sep 2013 11755/18797

ICA: Freeing Fourth Moments

• Create a matrix of fourth moment terms that would be

diagonal were the rows of H independent and diagonalize it

• A good candidate

– Good because it incorporates the energy in all rows of H – Where dij = E[ Sk hk

2 hi hj]

– i.e. D = E[hTh h hT]

• h are the columns of H
• Assuming h is real, else replace transposition with Hermition

50

         .. .. .. .. .. ..

23 22 21 13 12 11

d d d d d d D

23 Sep 2013 11755/18797

ICA: The D matrix

• Average above term across all columns of H

51

         .. .. .. .. .. ..

23 22 21 13 12 11

d d d d d d D

dij = E[ Sk hk

2 hi hj] = mj mi m k mk

h h h cols



2

) ( 1 H

jth component jth component Sum of squares

• f all components

Shk

2

hi hj Shk

2hi hj 23 Sep 2013 11755/18797

ICA: The D matrix

• If the hi terms were independent

– For i != j – Centered: E[hj] = 0  E[ Sk hk

2 hi hj]=0 for i != j

– For i = j

• Thus, if the hi terms were independent, dij = 0 if i != j
• i.e., if hi were independent, D would be a diagonal matrix

– Let us diagonalize D

52

         .. .. .. .. .. ..

23 22 21 13 12 11

d d d d d d D

dij = E[ Sk hk

2 hi hj] =

             

 

 

        

j k i k j i k i j j i k j i k

E E E E E E E E

, 2 3 3 2

h h h h h h h h h h

     

2 2 4 2

        

 

i k k i i k j i k

E E E E h h h h h h

mj mi m k mk

h h h cols



2

) ( 1 H

23 Sep 2013 11755/18797

Diagonalizing D

• Compose a fourth order matrix from X

– Recall: X = CM, H = BX = BCM

• B is what we’re trying to learn to make H independent

– Compose D’ = E[xT x x xT]

• Diagonalize D’ via Eigen decomposition

D’ = UUT

• B = UT

– That’s it!!!!

23 Sep 2013 11755/18797 53

B frees the fourth moment

D’ = UUT ; B = UT

• U is a unitary matrix, i.e. UTU = UUT = I (identity)
• H = BX = UTX
• h = UTx
• The fourth moment matrix of H is

E[hT h h hT] = E[xTUUTx UT x xTU] = E[xTx UT x xTU] = UT E[xTx xxT]U = UT D’ U = UT U  U T U = 

• The fourth moment matrix of H = UTX is Diagonal!!

54 23 Sep 2013 11755/18797

Overall Solution

• H = AM = BCM

– C is the (transpose of the) matrix of Eigen vectors of MMT

• X = CM
• A = BC = UTC

– B is the (transpose of the) matrix of Eigenvectors of X.diag(XTX).XT

55 23 Sep 2013 11755/18797

Independent Component Analysis

• Goal: to derive a matrix A such that the rows of AM are

independent

• Procedure:
• 1. “Center” M
• 2. Compute the autocorrelation matrix RMM of M
• 3. Compute whitening matrix C via Eigen decomposition

RMM = ESET, C = S-1/2ET

• 4. Compute X = CM
• 5. Compute the fourth moment matrix D’ = E[xTxxxT]
• 6. Diagonalize D’ via Eigen decomposition
• 7. D’ = UUT
• 8. Compute A = UTC
• The fourth moment matrix of H=AM is diagonal

– Note that the autocorrelation matrix of H will also be diagonal

56 23 Sep 2013 11755/18797

ICA by diagonalizing moment matrices

• The procedure just outlined, while fully functional, has

shortcomings

– Only a subset of fourth order moments are considered – There are many other ways of constructing fourth-order moment matrices that would ideally be diagonal

• Diagonalizing the particular fourth-order moment matrix we have chosen

is not guaranteed to diagonalize every other fourth-order moment matrix

• JADE: (Joint Approximate Diagonalization of Eigenmatrices),

J.F. Cardoso

– Jointly diagonalizes several fourth-order moment matrices – More effective than the procedure shown, but computationally more expensive

57 23 Sep 2013 11755/18797

Enforcing Independence

• Specifically ensure that the components of H are

independent

– H = AM

• Contrast function: A non-linear function that has a

minimum value when the output components are independent

• Define and minimize a contrast function

» F(AM)

• Contrast functions are often only approximations too..

23 Sep 2013 11755/18797 58

A note on pre-whitening

• The mixed signal is usually “prewhitened”

– Normalize variance along all directions – Eliminate second-order dependence

• Eigen decomposition MMT = ESET
• C = S-1/2ET
• Can use first K columns of E only if only K independent

sources are expected

– In microphone array setup – only K < M sources

• X = CM

– E[xixj] = E[xi]E[xj] = dij for centered signal

23 Sep 2013 11755/18797 59

The contrast function

• Contrast function: A non-linear function that

has a minimum value when the output components are independent

• An explicit contrast function
• With constraint : H = BX

– X is “whitened” M

23 Sep 2013 11755/18797 60

) ( ) ( ) ( h h H H H I

i i



 

Linear Functions

• h = Bx

– Individual columns of the H and X matrices – x is mixed signal, B is the unmixing matrix

23 Sep 2013 11755/18797 61

1

| | ) ( ) (

 

 B h B h

1 x h

P P



 x x x x d P P H ) ( log ) ( ) ( | | log ) ( ) ( B x h   H H

The contrast function

• Ignoring H(x) (Const)
• Minimize the above to obtain B

23 Sep 2013 11755/18797 62

) ( ) ( ) ( H h H H H I

i i



  | | log ) ( ) ( ) ( B x h H    H H I

i i

| | log ) ( ) ( B h H



 

i i

H J

An alternate approach

• Recall PCA
• M = WH, the columns of W must be

uncorrelated

• Leads to: minW||M –WTWM||2+trace(WWT)

– Error minimization framework to estimate W

• Can we arrive at an error minimization

framework for ICA

• Define an “Error” objective that represents

independence

23 Sep 2013 11755/18797 63

An alternate approach

• Definition of Independence – if x and y are

independent:

– E[f(x)g(y)] = E[f(x)]E[g(y)] – Must hold for every f() and g()!!

23 Sep 2013 11755/18797 64

An alternate approach

• Define g(H) = g(BX) (component-wise

function)

• Define f(H) = f(BX)

23 Sep 2013 11755/18797 65

g(h11) g(h12) . . . g(h21) g(h22) . . . . . . . . . . . . . . . f(h11) f(h12) f(h21) f(h22)

An alternate approach

• P = g(H) f(H)T = g(BX) f(BX)T

This is a square matrix

• Must ideally be
• Error = ||P-Q||F

2

23 Sep 2013 11755/18797 66

P11 P12 . . . P21 P22 . . . . . .





k jk ik ij

h f h g P ) ( ) (

. . .

P =

Q11 Q12 . . . Q21 Q22 . . .

Q =

j i h f h g Q

k l jl ik ij

 



) ( ) (





k il ik ii

h f h g Q ) ( ) (

An alternate approach

• Ideal value for Q
• If g() and h() are odd symmetric functions

Sjg(hij) = 0 for all i

– Since = Sjhij = 0 (H is centered) – Q is a Diagonal Matrix!!!

23 Sep 2013 11755/18797 67

. . . Q11 Q12 . . . Q21 Q22 . . .

Q =

j i h f h g Q

k l jl ik ij

 



) ( ) (





k il ik ii

h f h g Q ) ( ) (

An alternate approach

• Minimize Error
• Leads to trivial Widrow Hopf type iterative

rule:

23 Sep 2013 11755/18797 68

Diagonal  Q

T

g(BX)f(BX) P 

2

|| ||

F

error Q P  

T

g(BX)f(BX) E   Diag

T

EB B B   

Update Rules

• Multiple solutions under different

assumptions for g() and f()

• H = BX
• B = B +  DB
• Jutten Herraut : Online update

– DBij = f(hi)g(hj); -- actually assumed a recursive neural network

• Bell Sejnowski

– DB = ([BT]-1 – g(H)XT)

23 Sep 2013 11755/18797 69

Update Rules

• Multiple solutions under different

assumptions for g() and f()

• H = BX
• B = B +  DB
• Natural gradient -- f() = identity function

– DB = (I – g(H)HT)W

• Cichoki-Unbehaeven

– DB = (I – g(H)f(H)T)W

23 Sep 2013 11755/18797 70

What are G() and H()

• Must be odd symmetric functions
• Multiple functions proposed
• Audio signals in general

– DB = (I – HHT-Ktanh(H)HT)W

• Or simply

– DB = (I –Ktanh(H)HT)W

23 Sep 2013 11755/18797 71

      Gaussian sub is x ) tanh( Gaussian super is x ) tanh( ) ( x x x x x g

So how does it work?

• Example with instantaneous mixture of two

speakers

• Natural gradient update
• Works very well!

72 23 Sep 2013 11755/18797

Another example!

Input Mix Output

23 Sep 2013 11755/18797 73

Another Example

• Three instruments..

23 Sep 2013 11755/18797 74

The Notes

• Three instruments..

23 Sep 2013 11755/18797 75

ICA for data exploration

• The “bases” in PCA

represent the “building blocks”

– Ideally notes

• Very successfully used
• So can ICA be used to

do the same?

23 Sep 2013 11755/18797 76

ICA vs PCA bases

Non-Gaussian data ICA PCA

• Motivation for using ICA vs PCA
• PCA will indicate orthogonal directions
• f maximal variance
• May not align with the data!
• ICA finds directions that are

independent

• More likely to “align” with the data

23 Sep 2013 11755/18797 77

Finding useful transforms with ICA

• Audio preprocessing

example

• Take a lot of audio snippets

and concatenate them in a big matrix, do component analysis

• PCA results in the DCT bases
• ICA returns time/freq

localized sinusoids which is a better way to analyze sounds

• Ditto for images

– ICA returns localizes edge filters

23 Sep 2013 11755/18797 78

Example case: ICA-faces vs. Eigenfaces

ICA-faces Eigenfaces

23 Sep 2013 11755/18797 79

ICA for Signal Enhncement

• Very commonly used to enhance EEG signals
• EEG signals are frequently corrupted by

heartbeats and biorhythm signals

• ICA can be used to separate them out

23 Sep 2013 11755/18797 80

So how does that work?

• There are 12 notes in the segment, hence we

try to estimate 12 notes..

23 Sep 2013 11755/18797 81

PCA solution

• There are 12 notes in the segment, hence we

try to estimate 12 notes..

23 Sep 2013 11755/18797 82

So how does this work: ICA solution

• Better..

– But not much

• But the issues here?

23 Sep 2013 11755/18797 83

ICA Issues

• No sense of order

– Unlike PCA

• Get K independent directions, but does not have a notion
• f the “best” direction

– So the sources can come in any order – Permutation invariance

• Does not have sense of scaling

– Scaling the signal does not affect independence

• Outputs are scaled versions of desired signals in permuted
• rder

– In the best case – In worse case, output are not desired signals at all..

23 Sep 2013 11755/18797 84

What else went wrong?

• Notes are not independent

– Only one note plays at a time – If one note plays, other notes are not playing

• Will deal with these later in the course..

23 Sep 2013 11755/18797 85