Signal Decomposition via Distributed Optimization
Bennet Meyers and Stephen Boyd
Stanford University

SLAC AI Seminar, Menlo Park, April 23, 2021

Outline

Signal decomposition problem
Parameters and validation
Distributed solution method
PV data example2

Signals

I we consider vector valued signals that can include missing values

y1 , . . . , yT ∈ (R ∪ {?})p ,
y = y1 · · ·yTI partition indices into known and unknown values

K = {(t, i) | yt,i ∈ R},

U = {(t, i) | yt,i =?}I examples: financial data, energy production/load data, atmospheric andhydrospheric data

3123456780.60.4

kW

PV system number

Example: data from 8 PV systems

0.20.0time axis (minutes)
I 5 days of 1-minute measurements from 8 PV systems in the same geographic area
I white pixels denote missing values

4Decomposing a signal into components
I decomposition of signal y into K components
Ky = x1 + · · · + xK
K
I = means equal at known indexes t, i ∈ K
I x k ∈ RT ×p for k = 1, . . . , K are components with no missing data
I component x k comes from component class k
I example component classes– smooth– sparse– periodic– nonnegative– piecewise affine– Boolean with infrequent switching5

Brief comment on notation

I drop the k when refering to a general variable x ∈ RT ×p
I xt ∈ Rp for t = 1, . . . , T is a row vector
I xi ∈ RT for i = 1, . . . , p is a column vector
I xt,i ∈ R is a single entry

6

Estimating missing values

I missing values in y can be estimated from decomposition asUŷ = x 1 + · · · + x K
I basis of validation and parameter tuning method described later on

7Component classes
I component classes characterized by loss or implausibility functions

φk : RT ×p → R ∪ {∞},k = 1, . . . , K

I smaller φk (x) means more plausible x for class k
I infinite values encode constraints on components
I for statistical model of a component class, φ(x) is negative log-likelihood
I simple examples (with scale factor λ > 0):
λ– mean-square small class: φ(x) = Tp2t kxt k2

PPT −1
λ2– mean-square smooth class: φ(x) = (T −1)pt=1 kxt+1 − xt k2
0 xt,i ≥ 0 for all t, i
– nonnegative class: φ(x) =∞ otherwise8

Signal decomposition problem

I we choose decomposition to minimize total loss or implausibility
I signal decomposition (SD) problem:

minimize

φ1 (x 1 ) + · · · + φK (x K )Ksubject to y = x 1 + · · · + x K
I variables are components x 1 , . . . , x K
I we refer to a solution as an optimal signal decomposition

9

Solving the signal decomposition problem

I if all φk are convex, SD problem is convex, and so can be efficiently solved
I otherwise, we settle for an approximate solution

I our method is based on alternating directions method of multipliers (ADMM)– a distributed method that handles the component classes separately– easy to define new component classes– solves SD problem when it’s convex– approximately solves SD problem it’s not convex

10

ExampleK = 3 component classes:
I mean-square small, φ1 (x) = λ1Tp2t kxt k2
PT −1
λ2
2
I mean-square smooth, φ(x) =t=1 kxt+1 − xt k2(T −1)pP

I Boolean with infrequent switching,(φ(x) =P PT −1
λ3
it=1 |xt+1,i − xt,i |(T −1)p∞xt,i ∈ {0, 1} for all t, iotherwisei.e., rate of switching between values 0 and 1
I λi are positive weights; can take λ1 = 1

11Synthetic data

I scalar signal, i.e., y ∈ (R ∪ {∞})1×T
1 + x 2 + x 3 , then randomly make 10% of entries unknownI generate y as xtrue
true
true
1I entries of xtrue
are IID N (0, 0.152 )2
I entries of xtrue
are white noise passed through a low-pass filter
3I xtrueis realization of Markov chain on {0, 1} with probability 0.1 of transitions

0 → 1 or 1 → 0

12Synthetic data
component 1
0.25
0.00
−0.25

component 2
10−1

component 3
10signal, y20

missing data

−2025

50

75

100

12515017520013

Decomposition with λ2 = 25, λ3 = 0.5
component 1
trueestimate

0.25
0.00
−0.25

component 2true
estimate

0−1

component 3
1

true
estimate0

observed
missingtrue
estimatecomposed signal
200255075

100125150

17520014

Decomposition with λ2 = 500, λ3 = 5
component 1
1

true
estimate0

component 2
trueestimate

0−1

component 3
1

true
estimate0

observed
missingtrue
estimatecomposed signal
200255075

100125150

175

200

15

Outline

Signal decomposition problem
Parameters and validation
Distributed solution method
PV data example

16

Component class parameters
I class losses φk can have associated parameters, denoted φk (x k ; θk ), θk ∈ Θk
I some common examples– weight or scaling parameters φ(x; θ) = θ`(x), θ ∈ Θ = R++(often denoted with traditional symbol λ)T ×p– signal scaling parameters φ(x; θ) = φbool (x/θ) =⇒ x ∈ {0, θ}– constraint parameters φ(x; θ) = I(θ1 ≤ x ≤ θ2 )– basis parameters φ(x) = I(x = θα for some α ∈ Rd×n )I different parameters lead to different decompositions

17

Validating a decomposition
I randomly select test set T⊂ K and replace associated values in y with ?

I carry out decomposition using entries K \ TI decomposition yields estimates ŷt,i for (t, i) ∈ T
I quantify residuals or errors yt,i − ŷt,i , (t, i) ∈ T , with some metric,e.g. RMS or AA (average absolute)
I for more stable validation, process can repeat, e.g. k-fold cross validation orbootstrap sampling
I can be used to choose component classes and class parameters

18

Example
AA test error
I set λ1 = 1, search

0.45

1

10min error
best value within 2%

parameter space forbest λ2 and λ30.40

I randomly select 10%

0.35

λ3of the data point fortest set100

I decompose for each

0.30

parameter value(10 × 10 grid)I repeat 12 times andtake average error

0.25
10−1
100

101102

λ2
19

Final decomposition, λ2 = 21.5, λ3 = 0.278
component 1
trueestimate

0.25
0.00
−0.25

component 2
trueestimate

0−1

component 3
1

true
estimate0

composed signal

observedtrue
estimate

20025

50

75

100

125

150

175

20020

Outline

Signal decomposition problem
Parameters and validation
Distributed solution method
PV data example

21

Alternating direction direction of multipliers (ADMM)

I a method for solving convex optimization problems
I developed in 1970s, with roots in 1950s; modern treatment in Boyd et al. [2011]
I can be used as a heuristic for non-convex problems
I a distributed method, with different parts handled separately

22

SD via ADMM
I for iteration j = 1, . . .

(x k )j+1 := proxk (x k )j − u j ,ŷ

j+1:=KX

k = 1, . . . , K(x k )j+1k=1

2 j+1(ŷ− yt,i ), (t, i) ∈ KK t,iφk (x) + ρ2 kx − v k2F , proximal operator of φkj+1jut,i:= ut,i+I proxk (v ) = argminx

I ρ > 0 is an algorithm parameterjI ut,i are dual variables

23

Convergence and properties

I converges to (global) solution when all φk are convex, for any ρ > 0
I is a good heuristic in other cases, but choice of ρ can matter

I only need proximal operator for each component class
I first step can be carried out in parallel, for the k components
jI each component handled separately; coordination is via dual variables ut,i24

Proximal operator

ρ
proxφ (v ) = argmin φ(x) + kx − v k2F2x
I compromise between making φ(x) small and x near v
I when φ is an indicator function of a set C, proximal operator is projection onto C(and doesn’t depend on ρ)
I for many φ, proximal operator can be worked out analytically
I for others, can involve some computation

25

Examples
namemean-square smallφ(x)P
λ t kxt k22

average-absolute smallλP

mean-square small rth -order diff.

λP

non-negativeI(x ≥ 0)(v )+linear equality constraintI(Ax = b)

Boolean set

I(x ∈ {0, 1}T ×p )v − AT (AAT )−1 (Av − b)
(0 |vt,i | ≤ |vt,i − 1|1 otherwise

t kxt k1

2
i kDr xi k2

proxφ (v )ρ2λ+ρ vt,i vt,i − λ/ρ vt,i > λ/ρ
0|vt,i | < λ/ρvt,i + λ/ρ vt,i < −λ/ρ
−1TD
I + 2λDvr
ρ r

26

A less obvious example

I Boolean with infrequent switching
φ(x) =P PT −1
λ i t=1|xt+1,i − xt,i | xt,i ∈ {0, 1} for all t, i∞otherwise

I proximal operator can be evaluated by solving a graph shortest path problem

using dynamic programming, with cost O(T ) flops

27

Convergence example, convex casemin obj val

60

obj val50
403020100204060

80

100

iteration
28

Convergence example, non-convex case175

min obj val

obj val

150
125
100
7550
25020406080100

iteration29

Outline

Signal decomposition problem
Parameters and validation
Distributed solution method
PV data example

30

PV fleet outage detection
I recall data example of power signals from 8 PV systems (PV ‘fleet’) exposed tosimilar weather patterns
I want to automatically detect drops in system output that might be due to a

failure of a PV module or string of modules
I standard industry approach– make a physical model of each system– obtain local measurements of irradiance and temperature– compare actual to predicted for each system

I let’s try a purely data driven approach, using SD

31

The data set2345– 25% loss of power output insystem 6 during second-to-lastday– 50% loss of power output insystem 1 during final day

6systems

7

I artificially induce ‘failures’ in two

8I T=21600, p=99measurements from 9 systems

1fleet power data
I 15 days of 1-minutetime [minutes]
32

Data preprocessing

I scale each system data to about [0, 1] (use 95th percentile for UB, not max)
I take log10 and set zero values to ?
I taking log gives a multiplicative component model, instead of additive33

SD components
I residual: φ1 (x) = λ1 kxk2
F
I common clear sky component: smooth, equal across systems, daily periodic

PT −2– φ2 (x) = λ2 t=1 kxt − 2xt+1 + xt+2 k22– xi − xi+1 = 0, for i = 1, . . . , p − 1– xt − xt+1440 = 0, for t = 1, . . . , T − 1440
I common weather component: asymmetric distribution and equal across systems

P– φ3 (x) = λ3 t,i 1/2 |xt,i | + (τ − 1/2)xt,i– xt,i − xt,i+1 = 0, for all t and for i = 1, . . . , p − 1
I outage detector: non-positive, mostly zero, and mostly constant

P PT −1P– φ4 (x) = λ4 i t=1 |xt,i − xt+1,i | + λ5 t,i (−xt,i )– x 0– x1 = 0 (first row is the zero vector)
I T × p × K = 777, 636 variables to estimate!
34

Results: decomposition of system 1 (last 5 days)log power space
datacomp 1comp 20.0

power space
1

−2.5
0
1.5
1.0
0.5

0.0
−0.5
0.0

1

−2.50

comp 3comp 40−110
1.00

0.0

0.75

−0.2020004000

6000

02000

4000

6000
35

Results, outage detector component
Fleet loss analysis

PV system number

20.9

340.8

560.7

78

loss factor (unitless)

1.0

10.6

9time axis (minutes)
36

Naive approach, compare each to average
Fleet loss analysis

PV system number

21.0

340.8

560.6

78loss factor (unitless)

1.2

10.4

9time axis (minutes)
37

Software

I developing Python implementation:

https://github.com/bmeyers/optimal-signal-demixing/
I work in progress
I components defined as objects, proximal operators are attributes
I no requirement to understand ADMM or proximal operators to use!

38

Conclusions
signal decomposition via distributed optimization
I is interpretable
I provides a good way to describe prior knowledge about the signal
I is extensible
I is scalable to very large data sets
I does not require large training sets (or any labeled training data)

This material is based upon work supported by the U.S. Department ofEnergy’s Office of Energy Efficiency and Renewable Energy (EERE) under
the Solar Energy Technologies Office Award Number 34368.
39