jeffrey-hokanson / PSDR / 351

pending completion

Build # 351

Build Type

push

travis-ci-com

Committed by

jeffrey-hokanson

Commit Message

Fixing failing test for plot_envelope due to solver tolerances

Run Details

1205 of 1692 branches covered (71.22%)

Branch coverage included in aggregate %.

5 of 5 new or added lines in 1 file covered. (100.0%)

30 existing lines in 3 files now uncovered.

4513 of 5180 relevant lines covered (87.12%)

1.74 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

46.97

/psdr/opg.py

from __future__ import division, print_function

import numpy as np

from .subspace import SubspaceBasedDimensionReduction, ActiveSubspace

__all__ = ['OuterProductGradient']

def opg_grads(Z, fZ, kernel = None):
        if kernel is None:
                # Bandwidth from Xia 2007, [Li18, eq. 11.5] 
                bw = 2.34*len(Z)**(-1./(max(Z.shape[1], 3) +6))
                kernel = lambda dist: np.exp(-bw*dist**2/2.)
                 
        z_grads = np.zeros(Z.shape)                # estimated gradients in the transformed coordinates
        for i, zi in enumerate(Z):
                # This essentially fits the linear model by solving the normal equations
                A = np.zeros((Z.shape[1]+1, Z.shape[1]+1))
                b = np.zeros((Z.shape[1]+1))
                for j, zj in enumerate(Z):
                        h = np.hstack([1, zj - zi])
                        # TODO: Li mentions the kernel can be evaluated cheaper when Gaussian
                        kern = kernel(np.linalg.norm(zi - zj))
                        A += np.outer(h, h)*kern
                        b += h*fZ[j]*kern
                xx = np.linalg.solve(A, b)
                z_grads[i] = xx[1:]

        return z_grads


class OuterProductGradient(ActiveSubspace):
        r""" The Outer Product Gradient approach of Hardle and Stoker


        Given pairs of inputs :math:`\mathbf{x}_j` and outputs :math:`f(\mathbf{x}_j)`,
        this method estimates the active subspace by (1) estimating gradients from these samples
        and then (2) using the eigendecomposition of the outer product of these approximate gradients 
        The gradients are estimated by fitting a linear model at each input point :math:`\mathbf{x}_j`
        weighting the other points based on a kernel :math:`k` that decreases with increasing distance;
        i.e., at :math:`\mathbf{x}_j` we construct the linear model :math:`\mathbf{a}^\top \mathbf{x} + \beta`: 
        by solving a linear system:

        .. math::
                
                \min_{\mathbf{a}, \beta} 
                \left\|
                \begin{bmatrix}
                        \sqrt{k(\| \mathbf{x}_1 - \mathbf{x}_j\|)} & & & \\        
                        & \sqrt{k(\| \mathbf{x}_2 - \mathbf{x}_j\|)} & & \\
                        & & \ddots & \\
                        & & & \sqrt{k(\|\mathbf{x}_M - \mathbf{x}_j\|)}
                \end{bmatrix}
                \left[        
                \begin{bmatrix} 
                        \mathbf{x}_1^\top - \mathbf{x}_j^\top & 1 \\ 
                        \mathbf{x}_2^\top - \mathbf{x}_j^\top & 1 \\ 
                        \vdots & \vdots \\
                        \mathbf{x}_M^\top - \mathbf{x}_j^\top & 1 \\ 
                \end{bmatrix}
                \begin{bmatrix} \mathbf{a} \\ \beta \end{bmatrix}
                - 
                \begin{bmatrix}
                        f(\mathbf{x}_1) \\ 
                        f(\mathbf{x}_2) \\
                        \vdots \\
                        f(\mathbf{x}_M)
                \end{bmatrix}
                \right]
                \right\|_2.
                

        Our implementation largely follows [Li18]_.

        References
        ----------
        .. [HS89]  W. Hardle and T. M. Stoker.
                Investigating Smooth Multiple Regression by the Method of Average Derivatives
                J. Am. Stat. Ass. Vol 84, No 408, (1989). pp 986--995, DOI:0.1080/01621459.1989.10478863          
        .. [Li18] Bing Li.
                Sufficient Dimension Reduction: Methods and Applications with R.
                CRC Press, 2018.
        """
        def __init__(self, kernel = 'gaussian', standardize = True):
        
                # TODO: Bandwidth determination strategies
                assert kernel in ['gaussian',], "Invalid kernel provided"
                self.kernel = kernel
                self.standardize = standardize

        def __str__(self):
                return "<Outer Product Gradient>"

        def fit(self, X, fX):
                r""" 

                Parameters
                ----------
                X: array-like (M,m)
                        Location of function samples
                fX: array-like (M,)
                        Function samples where fX[j] = f(X[j])
                """

                X = np.atleast_2d(X)
                fX = np.atleast_1d(fX)

                if len(fX.shape) > 1:
                        fX = fX.flatten()
                        assert len(fX) == len(X), "Must provide one-dimensional inputs of same dimension as X"
                
                assert len(fX) == len(X), "Number of points X must match number of function samples fX"
                
                N = len(X)
                
                # Step 1: Standardize data
                if self.standardize:
                        mean = np.mean(X, axis = 0)
                        std = np.std(X, axis = 0)
                        Z = (X - mean)/std
                else:
                        Z = X

                # Step 2: Setup the kernel

                # Bandwidth from Xia 2007, [Li18, eq. 11.5] 
                bw = 2.34*len(X)**(-1./(max(X.shape[1], 3) +6))
                kernel = lambda dist: np.exp(-bw*dist**2/2.)

                # TODO: Implement both naive approach in Li as well as weighted LS
                # N.B.: weighted LS will be slower (but better conditioned) because it is solve a large LS
                # vs a small pos-definite linear system

                # Step 3: Estimate gradients
                z_grads = opg_grads(Z, fX, kernel)

                # Step 4: identify the active subspace
                ActiveSubspace.fit(self, z_grads)
                        
                if self.standardize:
                        # Due to the coordiante change introduced by standardization 
                        Dinv2 = np.diag(std**(-0.5))
                        U = Dinv2.dot(self.U)
                        Q, R = np.linalg.qr(U, mode = 'reduced')
                        self._U = Q
                        self._U = self._fix_subspace_signs_grads(self._U, z_grads)                
                


1	from __future__ import division, print_function	4✔
2
3	import numpy as np	4✔
4
5	from .subspace import SubspaceBasedDimensionReduction, ActiveSubspace	4✔
6
7	__all__ = ['OuterProductGradient']	4✔
8
9	def opg_grads(Z, fZ, kernel = None):	4✔
10	if kernel is None:	4!
11	# Bandwidth from Xia 2007, [Li18, eq. 11.5]
12	bw = 2.34len(Z)*(-1./(max(Z.shape[1], 3) +6))	4✔
13	kernel = lambda dist: np.exp(-bwdist*2/2.)	4✔
14
15	z_grads = np.zeros(Z.shape) # estimated gradients in the transformed coordinates	4✔
16	for i, zi in enumerate(Z):	4✔
17	# This essentially fits the linear model by solving the normal equations
18	A = np.zeros((Z.shape[1]+1, Z.shape[1]+1))	4✔
19	b = np.zeros((Z.shape[1]+1))	4✔
20	for j, zj in enumerate(Z):	4✔
21	h = np.hstack([1, zj - zi])	4✔
22	# TODO: Li mentions the kernel can be evaluated cheaper when Gaussian
23	kern = kernel(np.linalg.norm(zi - zj))	4✔
24	A += np.outer(h, h)*kern	4✔
25	b += hfZ[j]kern	4✔
26	xx = np.linalg.solve(A, b)	4✔
27	z_grads[i] = xx[1:]	4✔
28
29	return z_grads	4✔
30
31
32	class OuterProductGradient(ActiveSubspace):	4✔
33	r""" The Outer Product Gradient approach of Hardle and Stoker
34
35
36	Given pairs of inputs :math:`\mathbf{x}_j` and outputs :math:`f(\mathbf{x}_j)`,
37	this method estimates the active subspace by (1) estimating gradients from these samples
38	and then (2) using the eigendecomposition of the outer product of these approximate gradients
39	The gradients are estimated by fitting a linear model at each input point :math:`\mathbf{x}_j`
40	weighting the other points based on a kernel :math:`k` that decreases with increasing distance;
41	i.e., at :math:`\mathbf{x}_j` we construct the linear model :math:`\mathbf{a}^\top \mathbf{x} + \beta`:
42	by solving a linear system:
43
44	.. math::
45
46	\min_{\mathbf{a}, \beta}
47	\left\\|
48	\begin{bmatrix}
49	\sqrt{k(\\| \mathbf{x}_1 - \mathbf{x}_j\\|)} & & & \\
50	& \sqrt{k(\\| \mathbf{x}_2 - \mathbf{x}_j\\|)} & & \\
51	& & \ddots & \\
52	& & & \sqrt{k(\\|\mathbf{x}_M - \mathbf{x}_j\\|)}
53	\end{bmatrix}
54	\left[
55	\begin{bmatrix}
56	\mathbf{x}_1^\top - \mathbf{x}_j^\top & 1 \\
57	\mathbf{x}_2^\top - \mathbf{x}_j^\top & 1 \\
58	\vdots & \vdots \\
59	\mathbf{x}_M^\top - \mathbf{x}_j^\top & 1 \\
60	\end{bmatrix}
61	\begin{bmatrix} \mathbf{a} \\ \beta \end{bmatrix}
62	-
63	\begin{bmatrix}
64	f(\mathbf{x}_1) \\
65	f(\mathbf{x}_2) \\
66	\vdots \\
67	f(\mathbf{x}_M)
68	\end{bmatrix}
69	\right]
70	\right\\|_2.
71
72
73	Our implementation largely follows [Li18]_.
74
75	References
76	----------
77	.. [HS89] W. Hardle and T. M. Stoker.
78	Investigating Smooth Multiple Regression by the Method of Average Derivatives
79	J. Am. Stat. Ass. Vol 84, No 408, (1989). pp 986--995, DOI:0.1080/01621459.1989.10478863
80	.. [Li18] Bing Li.
81	Sufficient Dimension Reduction: Methods and Applications with R.
82	CRC Press, 2018.
83	"""
84	def __init__(self, kernel = 'gaussian', standardize = True):	4✔
85
86	# TODO: Bandwidth determination strategies
UNCOV 87	assert kernel in ['gaussian',], "Invalid kernel provided"	×
88	self.kernel = kernel	×
UNCOV 89	self.standardize = standardize	×
90
91	def __str__(self):	4✔
UNCOV 92	return "<Outer Product Gradient>"	×
93
94	def fit(self, X, fX):	4✔
95	r"""
96
97	Parameters
98	----------
99	X: array-like (M,m)
100	Location of function samples
101	fX: array-like (M,)
102	Function samples where fX[j] = f(X[j])
103	"""
104
UNCOV 105	X = np.atleast_2d(X)	×
UNCOV 106	fX = np.atleast_1d(fX)	×
107
UNCOV 108	if len(fX.shape) > 1:	×
UNCOV 109	fX = fX.flatten()	×
UNCOV 110	assert len(fX) == len(X), "Must provide one-dimensional inputs of same dimension as X"	×
111
112	assert len(fX) == len(X), "Number of points X must match number of function samples fX"	×
113
114	N = len(X)	×
115
116	# Step 1: Standardize data
117	if self.standardize:	×
UNCOV 118	mean = np.mean(X, axis = 0)	×
119	std = np.std(X, axis = 0)	×
120	Z = (X - mean)/std	×
121	else:
122	Z = X	×
123
124	# Step 2: Setup the kernel
125
126	# Bandwidth from Xia 2007, [Li18, eq. 11.5]
UNCOV 127	bw = 2.34len(X)*(-1./(max(X.shape[1], 3) +6))	×
128	kernel = lambda dist: np.exp(-bwdist*2/2.)	×
129
130	# TODO: Implement both naive approach in Li as well as weighted LS
131	# N.B.: weighted LS will be slower (but better conditioned) because it is solve a large LS
132	# vs a small pos-definite linear system
133
134	# Step 3: Estimate gradients
135	z_grads = opg_grads(Z, fX, kernel)	×
136
137	# Step 4: identify the active subspace
UNCOV 138	ActiveSubspace.fit(self, z_grads)	×
139
UNCOV 140	if self.standardize:	×
141	# Due to the coordiante change introduced by standardization
UNCOV 142	Dinv2 = np.diag(std**(-0.5))	×
UNCOV 143	U = Dinv2.dot(self.U)	×
UNCOV 144	Q, R = np.linalg.qr(U, mode = 'reduced')	×
UNCOV 145	self._U = Q	×
UNCOV 146	self._U = self._fix_subspace_signs_grads(self._U, z_grads)	×
147
148

jeffrey-hokanson / PSDR / 351

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous