**5. Nonlinear dimensionality reduction methods**

So far, we have investigated the use of linear dimensionality reduction methods (PCA and MDS) in extracting synergies. In this section we used nonlinear dimensionality reduction methods for the same purpose. The motivation to explore nonlinear methods was that physiologists who studied motor control have propounded that there were inherent nonlinearities in the human motor system. By using nonlinear methods we could probably achieve improved precision in reconstruction of natural movements. The nonlinear methods applied in this chapter are Isomap, local linear embedding (LLE), and kernel PCA (kPCA). The first two methods Isomap and LLE, are built on the framework of classical multidimensional scaling discussed in the previous section. kPCA is built on the framework of PCA.

### **5.1 Isomap**

Isomap is similar to PCA and MDS. Although Isomap does linear estimations in the data point neighborhoods, the synergies extracted are nonlinear because these small neighborhoods are stitched together without trying to maintain linearity.

**Angular**

**Samples Samples Samples**

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -1

**Synergy 4 Synergy 5 Synergy 6**

**Synergy 1 Synergy 2 Synergy 3**

Application of Linear and Nonlinear Dimensionality Reduction Methods 119




> -1 0 1


> -1 0 1


> -1 0 1


> 0 1







> 0 0.1 0.2


> 0 0.1 0.2

> 0 0.5

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -0.5

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -1



> -1 0 1

> -1 0 1

> -1 0 1

> -1 0 1

> -1 0 1

> -1 0 1

> -1 0 1

> > 0 1










> 0 0.5

**T\_MCP**

**IͲMCP IͲPIP**

**MͲMCP MͲPIP RͲMCP RͲPIP**

**PͲMCP PͲPIP**

**T\_MCP**

**IͲMCP IͲPIP**

**MͲMCP MͲPIP RͲMCP RͲPIP**

**PͲMCP PͲPIP**

**TͲIP**

**TͲIP**

**Samples Samples Samples**

Fig. 5. Six kinematic synergies obtained for subject 1 using Isomap. Each synergy is about 0.45 s in duration (39 samples at 86 Hz). Abbreviations: T, thumb; I, index finger; M, middle finger; R, ring finger; P, pinky finger; MCP, metacarpophalangeal joint; IP, interphalangeal

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -0.5

**(radians/sample)**

**Angular**

**velocities of ten joints**

**velocities of ten joints**

**(radians/sample)**



> -2 0 2


> -2 0 2


> -2 0 2


> -2 0 2





> -0.5 0 0.5


> -0.5 0 0.5


> -0.5 0 0.5

> -0.1 0 0.1

> > 0 0.5

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -1

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -0.5

joint; PIP, proximal IP joint.

The following were the steps involved in estimating nonlinear synergies using Isomap:


3. Find eigenvectors of *<sup>τ</sup>*(D), where *<sup>τ</sup>*(D) = <sup>−</sup>*HSH*/2, *Sij* = (*Dij*)<sup>2</sup> and *Hij* <sup>=</sup> *dij* <sup>−</sup> 1/*N*, where *N* is number of data points and *d* is Kronecker delta function.

In PCA we estimated the eigen values and eigen vectors of covariance of the data. Similarly, here, we took a nonlinear approach to preserve inter-point distances on the manifold. The matrix D is similar to covariance matrix in PCA. D can actually be thought of as the covariance matrix in higher dimensions. Since in an N-dimensional space, the dimensions are the data points, the covariance for a particular pair of dimensions is the distance between the data points that define those dimensions.

Although this method looks linear like PCA, the source of nonlinearity is the method in which inter-point distances are calculated. For Isomap, we do not use the Euclidean distances between the points. If we use, it becomes classical MDS discussed in previous section. Rather, we use those distances only for points considered neighbors. The rest of inter-point distances are calculated by finding the shortest path through the graph on the manifold using Floyd's algorithm ( Tenenbaum et al. (2000)). The goal of the Isomap is to preserve the geodesic distances rather than the euclidian distances. Geodesic distances are calculated by moving along the approximate nonlinear manifold with given data point and interpolation between them.

We used drtoolbox in MATLAB by van der Maaten et al. (2009) to perform Isomap on the angular-velocity matrix to extract synergies. Fig. 5showed the top six synergies extracted using this method. Similar to PCA, all the nonlinear methods also yield the nonlinear synergies in descending order of their significance. The synergies extracted using this method had more submovements when compared to those in PCA.

#### **5.2 Local Linear Embedding**

Locally Linear Embedding (LLE) as the name suggests, tries to find a nonlinear manifold by stitching together small linear neighborhoods ( Roweis & Saul (2000)). This is very similar to Isomap. The difference between the two algorithms is in how they do the stitching. Isomap does this by doing a graph traversal by preserving geodesic distances while LLE does it by finding a set of weights that perform local linear interpolations that closely approximate the data.

The following were the steps involved in estimating nonlinear synergies using LLE:


3. Given those weights, find new data points that minimize interpolation error in lower dimensional space

We used drtoolbox in MATLAB by van der Maaten et al. (2009) to perform LLE on the angular-velocity matrix to extract synergies. Fig. 6showed the top six synergies extracted using this method. Similar to Isomap, in this method also we found more submovements in synergies than those from PCA.

12 Will-be-set-by-IN-TECH

3. Find eigenvectors of *<sup>τ</sup>*(D), where *<sup>τ</sup>*(D) = <sup>−</sup>*HSH*/2, *Sij* = (*Dij*)<sup>2</sup> and *Hij* <sup>=</sup> *dij* <sup>−</sup> 1/*N*, where

In PCA we estimated the eigen values and eigen vectors of covariance of the data. Similarly, here, we took a nonlinear approach to preserve inter-point distances on the manifold. The matrix D is similar to covariance matrix in PCA. D can actually be thought of as the covariance matrix in higher dimensions. Since in an N-dimensional space, the dimensions are the data points, the covariance for a particular pair of dimensions is the distance between the data

Although this method looks linear like PCA, the source of nonlinearity is the method in which inter-point distances are calculated. For Isomap, we do not use the Euclidean distances between the points. If we use, it becomes classical MDS discussed in previous section. Rather, we use those distances only for points considered neighbors. The rest of inter-point distances are calculated by finding the shortest path through the graph on the manifold using Floyd's algorithm ( Tenenbaum et al. (2000)). The goal of the Isomap is to preserve the geodesic distances rather than the euclidian distances. Geodesic distances are calculated by moving along the approximate nonlinear manifold with given data point and interpolation between

We used drtoolbox in MATLAB by van der Maaten et al. (2009) to perform Isomap on the angular-velocity matrix to extract synergies. Fig. 5showed the top six synergies extracted using this method. Similar to PCA, all the nonlinear methods also yield the nonlinear synergies in descending order of their significance. The synergies extracted using this method

Locally Linear Embedding (LLE) as the name suggests, tries to find a nonlinear manifold by stitching together small linear neighborhoods ( Roweis & Saul (2000)). This is very similar to Isomap. The difference between the two algorithms is in how they do the stitching. Isomap does this by doing a graph traversal by preserving geodesic distances while LLE does it by finding a set of weights that perform local linear interpolations that closely approximate the

3. Given those weights, find new data points that minimize interpolation error in lower

We used drtoolbox in MATLAB by van der Maaten et al. (2009) to perform LLE on the angular-velocity matrix to extract synergies. Fig. 6showed the top six synergies extracted using this method. Similar to Isomap, in this method also we found more submovements in

The following were the steps involved in estimating nonlinear synergies using LLE:

2. Find weights that allow neighbors to interpolate original data accurately

The following were the steps involved in estimating nonlinear synergies using Isomap:

1. Define neighbors for each data point 2. Find D, a matrix of inter-point distances

points that define those dimensions.

**5.2 Local Linear Embedding**

1. Define neighbors for each data point

synergies than those from PCA.

them.

data.

dimensional space

*N* is number of data points and *d* is Kronecker delta function.

had more submovements when compared to those in PCA.

Fig. 5. Six kinematic synergies obtained for subject 1 using Isomap. Each synergy is about 0.45 s in duration (39 samples at 86 Hz). Abbreviations: T, thumb; I, index finger; M, middle finger; R, ring finger; P, pinky finger; MCP, metacarpophalangeal joint; IP, interphalangeal joint; PIP, proximal IP joint.

**5.3 Kernel PCA**

(2010a) for further details.

denoted **v**row:

Similarly, a synergy **s***<sup>j</sup>*

[*s j* <sup>1</sup>(1), ...,*s*

samples, then we obtain the following row vector:

and all their possible shifts with 1 ≤ *tjk* ≤ *T* − *ts*.

We add *T* − *ts* zeros after each *s*

with *tjk* zeros added before each *s*

*j*

*j i*

> *j i*

[0, ..., 0,*s j* <sup>1</sup>(1), ...,*s*

reconstruct the velocity profiles as in the following equation.

synergies were similar to those obtained from PCA.

**6. Reconstruction of natural and ASL movements**

Kernel PCA (kPCA) is an extension of PCA in a high-dimensional space ( Scholkopf et al. (1998)). A high-dimensional space is first constructed by using a kernel function. Instead of directly doing a PCA on the data, the kernel based high dimensional feature space is used as input. In this chapter, we have used a gaussian kernel function. Kernel PCA computes the principal eigenvectors of the kernel matrix, rather than those of the covariance matrix. A kernel matrix is similar to the inner product of the data points in the high dimensional space that is constructed using the kernel function. The application of PCA in the kernel space provides Kernel PCA the property of constructing nonlinear mappings. We used drtoolbox in MATLAB by van der Maaten et al. (2009) to perform kPCA on the angular-velocity matrix to extract synergies. Fig. 7showed the top six synergies extracted using this method. These

Application of Linear and Nonlinear Dimensionality Reduction Methods 121

The synergies extracted from linear and nonlinear dimensionality reduction methods were used in reconstruction of natural movements. *l*1-norm minimization was used to reconstruct natural and ASL movements from the extracted synergies. This method with illustrations was already presented in Vinjamuri et al. (2010a). We have included a brief explanation here for readability and for the sake of completeness. We ask the readers to refer Vinjamuri et al.

Briefly, these were the steps involved in *l*1-norm minimization algorithm that was used for reconstruction of natural and ASL movements. Let us assume for a subject *m* synergies were obtained. The duration of the synergies is *ts* samples (*ts* = 39 in this study). Consider an angular-velocity profile of the subject, {**v**(*t*),*t* = 1, ..., *T*}, where *T* (*T* = 82 in this study) represents the movement duration (in samples). This profile can be rewritten as a row vector,

**v**row = [*v*1(1), ..., *v*1(*T*), ..., *v*10(1), ..., *v*10(*T*)].

(·) can be rewritten as the following row vector:

length of the vector the same as that of **v**row. If the synergy is shifted in time by *tjk* (*tjk* ≤ *T* − *ts*)

Then we construct a matrix as shown in Fig. 8 consisting of the row vectors of the synergies

With the above notation, we are trying to achieve a linear combination of synergies that can

<sup>1</sup>(*ts*), 0, ..., 0, ...,

*j*

(1) and *T* − *ts* − *tjk* zeros added after each *s*

*j*

0, ..., 0,*s j* <sup>10</sup>(1), ...,*s*

*j* <sup>10</sup>(1), ...,*s*

*j*

<sup>10</sup>(*ts*), 0, ..., 0].

(*ts*) (*i* = 1, ..., 10) in the above vector in order to make the

<sup>10</sup>(*ts*), 0, ..., 0]

**v**row = **c***B* (8)

*j i* (*ts*).

<sup>1</sup>(*ts*), 0, ..., 0, ...,*s*

Fig. 6. Six kinematic synergies obtained for subject 1 using LLE. Each synergy is about 0.45 s in duration (39 samples at 86 Hz). Abbreviations: T, thumb; I, index finger; M, middle finger; R, ring finger; P, pinky finger; MCP, metacarpophalangeal joint; IP, interphalangeal joint; PIP, proximal IP joint.

#### **5.3 Kernel PCA**

14 Will-be-set-by-IN-TECH










> 0 0.1




> -0.2 0 0.2


> -0.1 0 0.1

> -0.1 0 0.1

> -0.1 0 0.1


> 0 0.2

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -0.1

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -0.2

**T\_MCP**

**IͲMCP IͲPIP**

**MͲMCP MͲPIP RͲMCP RͲPIP**

**PͲMCP PͲPIP**

**T\_MCP**

**IͲMCP IͲPIP**

**MͲMCP MͲPIP RͲMCP RͲPIP**

**PͲMCP PͲPIP**

**TͲIP**

**TͲIP**

**Synergy 1 Synergy 2 Synergy 3**



0.01 0.02 0.03

> -0.1 0 0.1

0 0.02 0.04

> -0.1 0 0.1

0 0.02 0.04

> -0.1 0 0.1

0.01 0.02 0.03

> 0 0.2





> -0.1 0 0.1

> -0.1 0 0.1

> -0.1 0 0.1


> -0.1 0 0.1

> > 0 0.1

**Samples Samples Samples**

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -0.2

**Synergy 4 Synergy 5 Synergy 6**

**Samples Samples Samples**

Fig. 6. Six kinematic synergies obtained for subject 1 using LLE. Each synergy is about 0.45 s in duration (39 samples at 86 Hz). Abbreviations: T, thumb; I, index finger; M, middle finger; R, ring finger; P, pinky finger; MCP, metacarpophalangeal joint; IP, interphalangeal joint; PIP,

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -0.1

**Angular**

**(radians/sample)**

**Angular**

**velocities of ten joints**

**velocities of ten joints**

**(radians/sample)**





> -0.1 0 0.1


> -0.1 0 0.1


> -0.1 0 0.1


> -0.1 0 0.1

> -0.1 0 0.1

> -0.2 0 0.2

> -0.1 0 0.1

> -0.2 0 0.2

> -0.1 0 0.1

> -0.2 0 0.2


> -0.2 0 0.2

0 0.05

proximal IP joint.

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -0.1

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -0.05

Kernel PCA (kPCA) is an extension of PCA in a high-dimensional space ( Scholkopf et al. (1998)). A high-dimensional space is first constructed by using a kernel function. Instead of directly doing a PCA on the data, the kernel based high dimensional feature space is used as input. In this chapter, we have used a gaussian kernel function. Kernel PCA computes the principal eigenvectors of the kernel matrix, rather than those of the covariance matrix. A kernel matrix is similar to the inner product of the data points in the high dimensional space that is constructed using the kernel function. The application of PCA in the kernel space provides Kernel PCA the property of constructing nonlinear mappings. We used drtoolbox in MATLAB by van der Maaten et al. (2009) to perform kPCA on the angular-velocity matrix to extract synergies. Fig. 7showed the top six synergies extracted using this method. These synergies were similar to those obtained from PCA.

#### **6. Reconstruction of natural and ASL movements**

The synergies extracted from linear and nonlinear dimensionality reduction methods were used in reconstruction of natural movements. *l*1-norm minimization was used to reconstruct natural and ASL movements from the extracted synergies. This method with illustrations was already presented in Vinjamuri et al. (2010a). We have included a brief explanation here for readability and for the sake of completeness. We ask the readers to refer Vinjamuri et al. (2010a) for further details.

Briefly, these were the steps involved in *l*1-norm minimization algorithm that was used for reconstruction of natural and ASL movements. Let us assume for a subject *m* synergies were obtained. The duration of the synergies is *ts* samples (*ts* = 39 in this study). Consider an angular-velocity profile of the subject, {**v**(*t*),*t* = 1, ..., *T*}, where *T* (*T* = 82 in this study) represents the movement duration (in samples). This profile can be rewritten as a row vector, denoted **v**row:

$$\mathbf{v}\_{\text{row}} = [v\_1(1), \dots, v\_1(T), \dots, v\_{10}(1), \dots, v\_{10}(T)].$$

Similarly, a synergy **s***<sup>j</sup>* (·) can be rewritten as the following row vector:

$$[s\_1^j(1), \ldots, s\_1^j(t\_s), 0, \ldots, 0, \ldots, s\_{10}^j(1), \ldots, s\_{10}^j(t\_s), 0, \ldots, 0].$$

We add *T* − *ts* zeros after each *s j i* (*ts*) (*i* = 1, ..., 10) in the above vector in order to make the length of the vector the same as that of **v**row. If the synergy is shifted in time by *tjk* (*tjk* ≤ *T* − *ts*) samples, then we obtain the following row vector:

$$\begin{aligned} [0, \ldots, 0, s\_1^{\dot{j}}(1), \ldots, s\_1^{\dot{j}}(t\_\mathbf{s}), 0, \ldots, 0, \ldots \\ 0, \ldots, 0, s\_{10}^{\dot{j}}(1), \ldots, s\_{10}^{\dot{j}}(t\_\mathbf{s}), 0, \ldots, 0] \end{aligned}$$

with *tjk* zeros added before each *s j i* (1) and *T* − *ts* − *tjk* zeros added after each *s j i* (*ts*).

Then we construct a matrix as shown in Fig. 8 consisting of the row vectors of the synergies and all their possible shifts with 1 ≤ *tjk* ≤ *T* − *ts*.

With the above notation, we are trying to achieve a linear combination of synergies that can reconstruct the velocity profiles as in the following equation.

$$\mathbf{v}\_{\text{row}} = \mathbf{c}B \tag{8}$$

…

… …

… …

… … … …

Fig. 8. (a) Synergies were rearranged as row vectors (b) A template matrix was formed by adding the time-shifted versions of synergies along the rows. Adapted from Vinjamuri et al.

[0, ..., *c*11, ..., 0, ..., *c*1*K*<sup>1</sup> , ..., 0, ..., *cm*1, ..., 0, ..., *cmKm* , ..., 0] with nonzero values *cjk* appearing at the (*T* − *ts* + 1)(*j* − 1) + *tjk*-th elements of **c**. The matrix *B* (shown in Fig. 8(b)) can be viewed as a bank or library of template functions with each row of *B* as a template. This bank can be overcomplete and contain linearly dependent subsets. Therefore, for a given movement profile **v**row and an overcomplete bank of template functions

We hypothesize that the strategy of central nervous system for dimensionality reduction in movement control is to use a small number of synergies and a small number of recruitments of these synergies for movement generation. Therefore, the coefficient vector **c** in (8) should be sparse, i.e., having a lot of zeros and only a small number of nonzero elements. Therefore,

The following was optimization problem that was used in selection of synergies in

where �·�<sup>2</sup> represents the *l*<sup>2</sup> norm or Euclidean norm of a vector and *λ* is a regulation

1

*<sup>λ</sup>* � **<sup>c</sup>***<sup>B</sup>* <sup>−</sup> **<sup>v</sup>**row �<sup>2</sup>

<sup>2</sup> (9)

…

Joint 1 Joint 2 Joint 10

… Joint 1 Joint 2 Joint 10

Application of Linear and Nonlinear Dimensionality Reduction Methods 123

…

(a)

(b)

Time-shifted versions of synergy 1

Time-shifted versions of synergy 2

*B*, there exists an infinite number of **c** satisfying (8).

reconstruction of a particular movement.

we seek the sparsest coefficient vector **c** such that **c***B* = **v**row.

Minimize � **c** �<sup>1</sup> +

(2010a).

parameter.

where **c** denotes

Synergy 1 Synergy 2

Fig. 7. Six kinematic synergies obtained for subject 1 using kPCA. Each synergy is about 0.45 s in duration (39 samples at 86 Hz). Abbreviations: T, thumb; I, index finger; M, middle finger; R, ring finger; P, pinky finger; MCP, metacarpophalangeal joint; IP, interphalangeal joint; PIP, proximal IP joint.

Fig. 8. (a) Synergies were rearranged as row vectors (b) A template matrix was formed by adding the time-shifted versions of synergies along the rows. Adapted from Vinjamuri et al. (2010a).

where **c** denotes

16 Will-be-set-by-IN-TECH


> -0.2 0 0.2

> -0.2 0 0.2

> -0.2 0 0.2

> -0.2 0 0.2

> -0.2 0 0.2

> -0.2 0 0.2

> -0.1 0 0.1

> -0.2 0 0.2

> > 0 0.2


> -0.2 0 0.2

> -0.1 0 0.1

> -0.1 0 0.1

> -0.1 0 0.1



> -0.1 0 0.1

> -0.1 0 0.1

> > 0 0.1

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -0.2

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -0.1

**T\_MCP**

**IͲMCP IͲPIP**

**MͲMCP MͲPIP RͲMCP RͲPIP**

**PͲMCP PͲPIP**

**T\_MCP**

**IͲMCP IͲPIP**

**MͲMCP MͲPIP RͲMCP RͲPIP**

**PͲMCP PͲPIP**

**TͲIP**

**TͲIP**

**Synergy 1 Synergy 2 Synergy 3**


> -0.1 0 0.1

> -0.2 0 0.2

> -0.5 0 0.5

> -0.5 0 0.5

> -0.5 0 0.5

> -0.2 0 0.2

> -0.2 0 0.2

> -0.2 0 0.2

> > 0 0.5

> > > -0.1 0 0.1


> -0.1 0 0.1

> -0.2 0 0.2

> -0.1 0 0.1

> -0.2 0 0.2

> -0.1 0 0.1

> -0.1 0 0.1

> -0.2 0 0.2

> > 0 0.2

**Samples Samples Samples**

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -0.5

**Synergy 4 Synergy 5 Synergy 6**

**Samples Samples Samples**

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -0.2

Fig. 7. Six kinematic synergies obtained for subject 1 using kPCA. Each synergy is about 0.45 s in duration (39 samples at 86 Hz). Abbreviations: T, thumb; I, index finger; M, middle finger; R, ring finger; P, pinky finger; MCP, metacarpophalangeal joint; IP, interphalangeal

**Angular**

**(radians/sample)**

**Angular**

**velocities of ten joints**

**velocities of ten joints**

**(radians/sample)**

0 0.1 0.2









> 0 1

0 0.05 0.1

> -0.2 0 0.2

> -0.1 0 0.1

> -0.2 0 0.2

> -0.1 0 0.1

> -0.1 0 0.1

> -0.1 0 0.1

> -0.1 0 0.1

> -0.1 0 0.1

> > 0 0.2

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -1

<sup>0</sup> <sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> -0.2

joint; PIP, proximal IP joint.

$$[0, \dots, c\_{11}, \dots, 0, \dots, c\_{1K\_1}, \dots, 0, \dots, c\_{m1}, \dots, 0, \dots, c\_{mK\_{m'}}, \dots, 0]$$

with nonzero values *cjk* appearing at the (*T* − *ts* + 1)(*j* − 1) + *tjk*-th elements of **c**. The matrix *B* (shown in Fig. 8(b)) can be viewed as a bank or library of template functions with each row of *B* as a template. This bank can be overcomplete and contain linearly dependent subsets. Therefore, for a given movement profile **v**row and an overcomplete bank of template functions *B*, there exists an infinite number of **c** satisfying (8).

We hypothesize that the strategy of central nervous system for dimensionality reduction in movement control is to use a small number of synergies and a small number of recruitments of these synergies for movement generation. Therefore, the coefficient vector **c** in (8) should be sparse, i.e., having a lot of zeros and only a small number of nonzero elements. Therefore, we seek the sparsest coefficient vector **c** such that **c***B* = **v**row.

The following was optimization problem that was used in selection of synergies in reconstruction of a particular movement.

$$\text{Minimize} \quad \|\mathbf{c}\|\_1 + \frac{1}{\lambda} \left\| \|\mathbf{c}\mathbf{B} - \mathbf{v}\_{\text{row}}\|\_2^2 \tag{9}$$

where �·�<sup>2</sup> represents the *l*<sup>2</sup> norm or Euclidean norm of a vector and *λ* is a regulation parameter.

these synergies to reconstruct natural movements that were similar to activities of daily living. To broaden the applicability of synergies we also tested them on ASL postural movements which are different movements when compared to natural reach and grasp movements. We employed four types of dimensionality reduction methods: (1) PCA and MDS (2) Isomap (3) LLE (4) kPCA. PCA is a well known linear dimensionality reduction method. Two widely used PCA implementations (Covariance and SVD) were presented and relevant MATLAB codes were provided. Classical MDS is very similar to PCA but operates on a distance matrix. This was introduced as Isomap and LLE both work on a similar framework. Isomap and LLE are both neighborhood graph based nonlinear dimensionality reduction methods. The difference between Isomap and LLE is that the former is a global dimensionality reduction technique where as the former as the name suggests is a local linear technique. kPCA is similar to Isomap as it is a global dimensionality reduction method, but the uniqueness of the method is in using kernel tricks to transform the input data to higher dimensional space.

Application of Linear and Nonlinear Dimensionality Reduction Methods 125

Subject PCA Isomap LLE kPCA

Table 1. Mean reconstruction errors (± standard deviation) in natural movements

Subject PCA Isomap LLE kPCA

Table 2. Mean reconstruction errors (± standard deviation) in ASL postural movements

employed in this chapter.

Thus the reader was given the opportunity to sample different varieties of dimensionality reduction methods. Quantitative and qualitative comparison of the results obtained from reconstruction follows but the verdict is that PCA outperformed the nonlinear methods

The results from the reconstructions reveal that nonlinear techniques do not outperform the traditional PCA for both natural movements as well as ASL postural movements. The

1 0.2449 ± 0.0986 0.3480 ± 0.0988 0.3867±0.1052 0.3654 ± 0.0865 2 0.2782 ± 0.0915 0.3280 ± 0.0891 0.3661±0.0963 0.3291 ± 0.0745 3 0.2677 ± 0.0866 0.3110 ± 0.0764 0.3527±0.0881 0.3321 ± 0.0666 4 0.2321 ± 0.0645 0.3450 ± 0.0662 0.3667±0.1123 0.3540 ± 0.0831 5 0.2899 ± 0.0713 0.3780 ± 0.0843 0.3889±0.0792 0.3624 ± 0.0934 6 0.2569 ± 0.0456 0.3200 ± 0.0894 0.3457±0.0856 0.3654 ± 0.0678 7 0.2319 ± 0.0521 0.3890 ± 0.0663 0.3999±0.0996 0.3889 ± 0.0610 8 0.2789 ± 0.0486 0.3671 ± 0.0713 0.3897±0.1152 0.3764 ± 0.0885 9 0.2999 ± 0.0816 0.3335 ± 0.0465 0.3867±0.1672 0.3478 ± 0.0585 10 0.2569 ± 0.0936 0.3146 ± 0.0548 0.3547±0.0082 0.3354 ± 0.0635

1 0.0797 ± 0.0596 0.3240 ± 0.0769 0.3544±0.0954 0.2954 ± 0.0851 2 0.0845 ± 0.0781 0.3780 ± 0.0519 0.3998±0.0765 0.3444 ± 0.0658 3 0.0678 ± 0.0396 0.3180 ± 0.0554 0.3234±0.0834 0.2654 ± 0.0767 4 0.0590 ± 0.0432 0.3680 ± 0.0663 0.3879±0.0645 0.3214 ± 0.0873 5 0.0997 ± 0.0831 0.2640 ± 0.0734 0.2950±0.0790 0.2230 ± 0.0667 6 0.0661 ± 0.0976 0.3240 ± 0.0545 0.3754±0.0531 0.2875 ± 0.0931 7 0.0598 ± 0.0542 0.4140 ± 0.0787 0.4344±0.0632 0.3895 ± 0.0696 8 0.0732 ± 0.0348 0.3540 ± 0.0989 0.3954±0.0854 0.3123 ± 0.0555 9 0.0814 ± 0.0212 0.2290 ± 0.0823 0.3100±0.0991 0.2074 ± 0.0651 10 0.0883 ± 0.0443 0.2490 ± 0.0665 0.2780±0.0799 0.2189 ± 0.0799

Using the above optimization algorithm, the synergies extracted from four methods (PCA, Isomap, LLE, and kPCA) were used in reconstruction of natural movements and ASL postural movements. The reconstruction errors were calculated using the methods in Vinjamuri et al. (2010a). Figures 9 and 10 showed the comparison between the four dimensionality reduction methods. Fig. 9 showed the reconstruction errors for 100 natural movements and Fig. 10 showed the reconstruction errors for 36 ASL postural movements for all four methods. It is observed that PCA still has the best overall performance when compared with the novel nonlinear methods. Fig. 11 showed one of the best reconstructions by all four methods. Tables 1 and 2 summarize the reconstruction results obtained for all ten subjects.

Fig. 9. The reconstruction errors for 100 natural movements with four dimensionality reduction methods PCA, Isomap, LLE, and kPCA. PCA outperformed other three nonlinear methods.

Fig. 10. The reconstruction errors for 36 ASL postural movements with four dimensionality reduction methods PCA, Isomap, LLE, and kPCA. All methods performed poorly in reconstruction of ASL movements when compared to natural movements. PCA performed better than other three nonlinear methods.

#### **7. Summary**

In this chapter we applied linear and nonlinear dimensionality reduction methods to extract movement primitives or synergies from rapid reach and grasp movements. We then used 18 Will-be-set-by-IN-TECH

Using the above optimization algorithm, the synergies extracted from four methods (PCA, Isomap, LLE, and kPCA) were used in reconstruction of natural movements and ASL postural movements. The reconstruction errors were calculated using the methods in Vinjamuri et al. (2010a). Figures 9 and 10 showed the comparison between the four dimensionality reduction methods. Fig. 9 showed the reconstruction errors for 100 natural movements and Fig. 10 showed the reconstruction errors for 36 ASL postural movements for all four methods. It is observed that PCA still has the best overall performance when compared with the novel nonlinear methods. Fig. 11 showed one of the best reconstructions by all four methods. Tables

10 20 30 40 50 60 70 80 90 100

5 10 15 20 25 30 35

Fig. 10. The reconstruction errors for 36 ASL postural movements with four dimensionality reduction methods PCA, Isomap, LLE, and kPCA. All methods performed poorly in reconstruction of ASL movements when compared to natural movements. PCA performed

In this chapter we applied linear and nonlinear dimensionality reduction methods to extract movement primitives or synergies from rapid reach and grasp movements. We then used

Fig. 9. The reconstruction errors for 100 natural movements with four dimensionality reduction methods PCA, Isomap, LLE, and kPCA. PCA outperformed other three nonlinear

1 and 2 summarize the reconstruction results obtained for all ten subjects.

0 0.5 1

PCA

LLE

kPCA

PCA

0 0.5 1

0 0.5 1

0 0.5 1

0 0.5 1

better than other three nonlinear methods.

LLE

kPCA

**7. Summary**

 Isomap

methods.

 Isomap

0 0.5 1

0 0.5 1

0 0.5 1 these synergies to reconstruct natural movements that were similar to activities of daily living. To broaden the applicability of synergies we also tested them on ASL postural movements which are different movements when compared to natural reach and grasp movements. We employed four types of dimensionality reduction methods: (1) PCA and MDS (2) Isomap (3) LLE (4) kPCA. PCA is a well known linear dimensionality reduction method. Two widely used PCA implementations (Covariance and SVD) were presented and relevant MATLAB codes were provided. Classical MDS is very similar to PCA but operates on a distance matrix. This was introduced as Isomap and LLE both work on a similar framework. Isomap and LLE are both neighborhood graph based nonlinear dimensionality reduction methods. The difference between Isomap and LLE is that the former is a global dimensionality reduction technique where as the former as the name suggests is a local linear technique. kPCA is similar to Isomap as it is a global dimensionality reduction method, but the uniqueness of the method is in using kernel tricks to transform the input data to higher dimensional space.


Table 1. Mean reconstruction errors (± standard deviation) in natural movements



Thus the reader was given the opportunity to sample different varieties of dimensionality reduction methods. Quantitative and qualitative comparison of the results obtained from reconstruction follows but the verdict is that PCA outperformed the nonlinear methods employed in this chapter.

The results from the reconstructions reveal that nonlinear techniques do not outperform the traditional PCA for both natural movements as well as ASL postural movements. The

reconstruction errors were more for ASL postural movements when compared to those of natural movements for all methods. The reconstruction errors were in general larger for Isomap and LLE when compared with PCA and kPCA and of course PCA had outstanding performance for more than 90% of the tasks. van der Maaten et al. ( van der Maaten et al. (2009)) also found that nonlinear methods performed well on specific data sets but could not perform better than PCA for real world tasks. For example, for the Swiss roll data set that contains points that lie on a spiral like two dimensional manifold within a three dimensional space, several nonlinear techniques such as Isomap, LLE were able to find the two dimensional planar embedding, but linear techniques like PCA failed to find so. The reasons for two nonlinear methods Isomap and LLE to perform poorly in this study, might be that they relied neighborhood graphs. Moreover LLE might have been biased to local properties that do not necessarily follow the global properties of high dimensional data. It was surprising to see that Isomap, being a global dimensionality reduction technique performed poorly when compared to LLE for natural movements. kPCA performed better than Isomap and LLE, but kPCA does suffer from the limitation of selection of ideal kernel. The selection of gaussian kernel in this study might not be favorable in extracting the kinematic synergies in this study. In conclusion, although there are numerical advantages and disadvantages with both linear and nonlinear dimensionality reduction methods, PCA seemed to generalize and

Application of Linear and Nonlinear Dimensionality Reduction Methods 127

This work was supported by the NSF grant CMMI-0953449, NIDRR grant H133F100001. Special thanks to Laurens van der Maaten for guidance with the dimensionality reduction toolbox, and Prof. Dan Ventura (Brigham Young University) for helpful notes on comparison of LLE and Isomap. Thanks to Stephen Foldes for his suggestions with formatting. Thanks to Mr. Oliver Kurelic for his guidance and help through the preparation of the manuscript.

Bernabucci, I., Conforto, S., Capozza, M., Accornero, N., Schmid, M. & D'Alessio, T. (2007). A

Bernstein, N. (1967). *The Co-ordination and Regulation of Movements*, Pergamon Press, Oxford,

Braido, P. & Zhang, X. (2004). Quantitative analysis of finger motion coordination in hand

Cole, K. J. & Abbs, J. H. (1986). Coordination of three-joint digit movements for rapid

d'Avella, A., Portone, A., Fernandez, L. & Lacquaniti, F. (2006). Control of

d'Avella, A., Saltiel, P. & Bizzi, E. (2003). Combinations of muscle synergies in the contruction

Flash, T. & Hochner, B. (2005). Motor primitives in vertebrates and invertebrates, *Current*

Grinyagin, I. V., Biryukova, E. V. & Maier, M. A. (2005). Kinematic and dynamic synergies of human precision-grip movements, *Journal of Neurophysiology* 94(4): 2284–2294.

manipulative and gestic acts, *Human Movement Science* 22: 661–678.

finger-thumb grasp, *Journal of Neurophysiology* 55(6): 1407–1423.

of a natural motor behavior, *Nature Neuroscience* 6(3): 300–308.

*of NeuroEngineering and Rehabilitation* 4(33): 1–17.

*Opinion in Neurobiology* 15(6): 660–666.

biologically inspired neural network controller for ballistic arm movements, *Journal*

fast-reaching movements by muscle synergy combinations, *Journal of Neuroscience*

perform well on the real world data.

**8. Acknowledgments**

**9. References**

UK.

26(30): 7791–7810.

Fig. 11. An example reconstruction (in black) of a natural movement (in red) for task 24 when subject 1 was grasping an object.

reconstruction errors were more for ASL postural movements when compared to those of natural movements for all methods. The reconstruction errors were in general larger for Isomap and LLE when compared with PCA and kPCA and of course PCA had outstanding performance for more than 90% of the tasks. van der Maaten et al. ( van der Maaten et al. (2009)) also found that nonlinear methods performed well on specific data sets but could not perform better than PCA for real world tasks. For example, for the Swiss roll data set that contains points that lie on a spiral like two dimensional manifold within a three dimensional space, several nonlinear techniques such as Isomap, LLE were able to find the two dimensional planar embedding, but linear techniques like PCA failed to find so. The reasons for two nonlinear methods Isomap and LLE to perform poorly in this study, might be that they relied neighborhood graphs. Moreover LLE might have been biased to local properties that do not necessarily follow the global properties of high dimensional data. It was surprising to see that Isomap, being a global dimensionality reduction technique performed poorly when compared to LLE for natural movements. kPCA performed better than Isomap and LLE, but kPCA does suffer from the limitation of selection of ideal kernel. The selection of gaussian kernel in this study might not be favorable in extracting the kinematic synergies in this study. In conclusion, although there are numerical advantages and disadvantages with both linear and nonlinear dimensionality reduction methods, PCA seemed to generalize and perform well on the real world data.

#### **8. Acknowledgments**

20 Will-be-set-by-IN-TECH

0 0.02

0 0.05

0 0.05

0 0.05

0 0.05

0 0.05

0 0.05

0 0.05

0 0.05

> 0 0.1

0 0.01

0 0.05

0 0.05

0 0.05

0 0.05

0 0.05

0 0.05

0 0.05

0 0.05

> 0 0.1

Fig. 11. An example reconstruction (in black) of a natural movement (in red) for task 24 when

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.02

**Isomap**

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.1

**Samples**

**kPCA**

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.01

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.1

**Samples**

**T\_MCP**

**T\_MCP**

**IͲMCP**

**IͲPIP**

**MͲMCP**

**MͲPIP RͲMCP**

**RͲPIP**

**PͲMCP**

**PͲPIP**

**TͲIP**

**IͲMCP**

**IͲPIP**

**MͲMCP**

**MͲPIP RͲMCP**

**RͲPIP**

**PͲMCP**

**PͲPIP**

**TͲIP**

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.01

**PCA**

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.1

**Samples**

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.02

**LLE**

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.05

<sup>10</sup> <sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> -0.1

**Samples**

0 0.01

0 0.05

0 0.05

0 0.05

0 0.05

0 0.05

0 0.05

0 0.05

**velocities of ten**

**joints**

**(radians/sample)**

0 0.05

0 0.1 **Angular**

> 0 0.02

> 0 0.05

> 0 0.05

> 0 0.05

> 0 0.05

> 0 0.05

> 0 0.05

> 0 0.05

> 0 0.05

> > 0 0.1

**Angular**

subject 1 was grasping an object.

**velocities of ten**

**joints**

**(radians/sample)**

This work was supported by the NSF grant CMMI-0953449, NIDRR grant H133F100001. Special thanks to Laurens van der Maaten for guidance with the dimensionality reduction toolbox, and Prof. Dan Ventura (Brigham Young University) for helpful notes on comparison of LLE and Isomap. Thanks to Stephen Foldes for his suggestions with formatting. Thanks to Mr. Oliver Kurelic for his guidance and help through the preparation of the manuscript.

#### **9. References**


**0**

**7**

**Components Analysis**

<sup>1</sup>*Okayama University of Science*

<sup>2</sup>*Okayama University*

*Japan*

**Acceleration of Convergence of the Alternating**

**Least Squares Algorithm for Nonlinear Principal**

Masahiro Kuroda1, Yuichi Mori1, Masaya Iizuka2 and Michio Sakakihara1

Principal components analysis (PCA) is a popular descriptive multivariate method for handling quantitative data. In PCA of a mixture of quantitative and qualitative data, it requires quantification of qualitative data to obtain optimal scaling data and use ordinary PCA. The extended PCA including such quantification is called *nonlinear PCA*, see Gifi [Gifi, 1990]. The existing algorithms for nonlinear PCA are PRINCIPALS of Young et al. [Young et al., 1978] and PRINCALS of Gifi [Gifi, 1990] in which the alternating least squares (ALS) algorithm is utilized. The algorithm alternates between quantification of qualitative data and

In the application of nonlinear PCA for very large data sets and variable selection problems, many iterations and much computation time may be required for convergence of the ALS algorithm, because its speed of convergence is linear. Kuroda et al. [Kuroda et al., 2011] proposed an acceleration algorithm for speeding up the convergence of the ALS algorithm using the vector *ε* (v*ε*) algorithm of Wynn [Wynn, 1962]. During iterations of the v*ε* accelerated ALS algorithm, the v*ε* algorithm generates an accelerated sequence of optimal scaling data estimated by the ALS algorithm. Then the v*ε* accelerated sequence converges faster than the original sequence of the estimated optimal scaling data. In this paper, we use PRINCIPALS as the ALS algorithm for nonlinear PCA and provide the v*ε* acceleration for PRINCIPALS (v*ε*-PRINCIPALS). The computation steps of PRINCALS are given in Appendix A. As shown

in Kuroda et al. [Kuroda et al., 2011], the v*ε* acceleration is applicable to PRINCALS.

v*ε*-PRINCIPALS. In Section 6, we present our concluding remarks.

The paper is organized as follows. We briefly describe nonlinear PCA of a mixture of quantitative and qualitative data in Section 2, and describe PRINCIPALS for finding least squares estimates of the model and optimal scaling parameters in Section 3. Section 4 presents the procedure of v*ε*-PRINCIPALS that adds the v*ε* algorithm to PRINCIPALS for speeding up convergence and demonstrate the performance of the v*ε* acceleration using numerical experiments. In Section 5, we apply v*ε*-PRINCIPALS to variable selection in nonlinear PCA. Then we utilize modified PCA (M.PCA) approach of Tanaka and Mori [Tanaka and Mori, 1997] for variable selection problems and give the variable selection procedures in M.PCA of qualitative data. Numerical experiments examine the the performance and properties of

computation of ordinary PCA of optimal scaling data.

**1. Introduction**

