

Abstract-In this paper, a direct preconditioner is proposed for coupling matrix reconfiguration of bandpass filters with irregular couplings using the coefficient-parameter continuation method. With the preconditioner, the continuation method can be effectively used for filter synthesis and model extraction of bandpass filters with irregular couplings. In dealing with practical filters that contain a coupling topology with regular cross couplings introduced intentionally or unintentionally, the proposed preconditioner, being used as the start system in the continuation method, can swiftly lead to the solution that matches the physical realization. The preconditioner is physics-based and can be easily obtained by reconfiguring the regular coupling topology with the conventional method. The continuation process with a coefficient-parameter homotopy for coupling matrix reconfiguration is formulated for the first time, which transforms a coupling matrix in canonical form to that in the desired one with irregular couplings. The convergence of the solution and a bootstrapping approach for dealing with multiple irregular couplings are also discussed. The method provides an enabling tool for robot automatic tuning (RAT) of advanced microwave bandpass filters. Three practical examples are presented, including one synthesis example and two model extraction examples for RAT, demonstrating the effectiveness of the reconfiguration method for practical applications.

Index Terms-Coupling matrix, filter reconfiguration, homotopy continuation, irregular topologies, microwave filter tuning, preconditioner.

I. INTRODUCTION

UNING of high-performance microwave filters is inevitable due to stringent electrical specifications. The labor-intensive manual tuning process, being time-consuming and costly, motivates the recent development of computer-aided tuning (CAT) and robot automatic tuning (RAT), for which coupling matrix [START_REF] Cameron | General coupling matrix synthesis methods for Chebyshev filtering functions[END_REF] is used as the bridge between the physical model and the circuit model of the filter realization in the filter domain. In a robust RAT system, extraction of a legitimate circuit model contains three major Manuscript received May 7, 2020. The work described in this paper was substantially supported by the Hong Kong Ph.D Fellowship. (Corresponding author: Ke-Li Wu.)

Yan Zhang and Ke-Li Wu are with the Department of Electronic Engineering, The Chinese University of Hong Kong, Shatin, Hong Kong (e-mail: yzhang@link.cuhk.edu.hk; klwu@cuhk.edu.hk). Fabien Seyfert is with INRIA, Sophia-Antipolis, France, (e-mail: fabien.seyfert@inria.fr). steps: 1) to obtain the filtering function characterized by three polynomials from the measured response; 2) to transform the filtering function into the circuit model (coupling matrix) in a canonical form; and 3) to reconfigure the coupling matrix in the canonical form to that in the coupling topology that best fits the physical realization. The third step is usually conducted through a sequence of orthogonal transformations. The subsequent tuning strategy of the filter will depend on the difference between the extracted and the target coupling matrices. Among the three steps, the system identification problem in the first step can be addressed using well-developed methods, such as vector fitting [START_REF] Hu | A generalized coupling matrix extraction technique for bandpass filters with uneven-Qs[END_REF] and the recently proposed model-based vector fitting (MVF) methods [START_REF] Zhao | Model-based vector-fitting method for circuit model extraction of coupled-resonator diplexers[END_REF]. As long as a set of sensible filtering functions are obtained, a unique coupling matrix in a canonical form can be derived straightforwardly in the second step. However, coupling matrix reconfiguration in the third step using well-established orthogonal transformations [START_REF] Cameron | Microwave Filters for Communication Systems: Fundamentals Design and Applications[END_REF] is limited to the conventional regular coupling topologies, such as those constructed by cascaded trisections (CTs) and cascaded quartets (CQs) [START_REF] Tamiazzo | An analytical technique for the synthesis of cascaded N-tuplets cross-coupled resonators microwave filters using matrix rotations[END_REF].

The demands for dealing with bandpass filters in irregular coupling topologies are frequently encountered in many practical scenarios not only for model extraction but also filter synthesis. Such filters include those with 1) unintentional parasitic cross couplings; 2) intentional but unconventional cross couplings for enhancing the rejection; 3) coupling topologies for evenly distributing power among the internal nodes for high power applications [START_REF] Senior | Optimum network topologies for high power microwave filters[END_REF]; and 4) a weak dispersion effect, which can be described by parasitic cross couplings.

To reconfigure the coupling matrix of a bandpass filter with an irregular coupling topology, several algorithms are available in the literature. In [START_REF] Cameron | An exhaustive approach to the coupling matrix synthesis problem application to the design of high degree asymmetric filters[END_REF], an exhaustive method is applied to solve a set of polynomial equations by searching for an orthogonal matrix that enforces the non-existing couplings to be zero after similarity transformation. An analytic elimination can be imposed on a set of multiple variable polynomial system to obtain a Gröbner basis for the target coupling topology, with which the orthogonal transformation matrix for obtaining the coupling matrix of an irregular coupling topology can be obtained straightforwardly. Although multiple possible solutions can be found with the Gröbner basis method [START_REF] Seyfert | Dedale-HF[END_REF], it is difficult for a RAT application to stick to the target solution that is pertinent to the physical realization.

Some numerical reconfiguration methods using non-linear

T

optimization have also been developed. For example, Gauss-Newton method is used in [START_REF] Macchiarella | A powerful tool for the synthesis of prototype filters with arbitrary topology[END_REF] to solve a nonlinear least square problem. A new cost function based on the reflection and transmission coefficients at system zeros and poles of a trial coupling matrix is proposed in [START_REF] Atia | Synthesis of general topology multiple coupled resonator filters by optimization[END_REF] and [START_REF] Amari | Synthesis of cross-coupled resonator filters using an analytical gradient-based optimization technique[END_REF], which is minimized to find the coupling matrix with the desired coupling topology. Although these methods have shown some encouraging results for filter synthesis problems, the long computing time and poor convergence to the desired solution make them unfavorable for RAT applications.

To obtain global convergence, the genetic algorithm that is incorporated with local gradient-based optimization is proposed in [START_REF] Uhm | Synthesis of resonator filters with arbitrary topology using hybrid method[END_REF]. The cost function is based on the difference of the polynomial coefficients of the desired and the derived ones from a trial coupling matrix in the desired topology. Multiple possible solutions can be found using randomly generated initial populations in multiple runs. Very recently, an isospectral flow method that avoids direct optimization of the entries in a coupling matrix is proposed [START_REF] Pflüger | Coupling matrix extraction and reconfiguration using a generalized isospectral flow method[END_REF]. The method gradually adjusts the basis of the vector space pointing to the direction that better fulfills the desired coupling matrix structure. To deal with the multiple solution issue, multiple runs of iterations with random initial values have to be applied exhaustively.

The aforementioned methods have to face the following prohibitive issues for practical use in a RAT process: 1) slow convergence for high order filters; 2) easy to be trapped to a local minimum; and awkwardly, 3) likely falling into an undesired solution. Having said that, this work is the first attempt to overcome these predicaments in the reconfiguration of coupling matrix with an irregular coupling topology.

For RAT applications, it is desirable to extract the coupling matrix from the measured data that 1) matches the physical coupling topology of the filter; 2) is the right solution aligned with the target coupling matrix among multiple possible solutions; and 3) is acquired swiftly to accommodate the real-time tuning. To the best of the authors' knowledge, none of the existing reconfiguration approaches is tailored for these attributes. It is particularly true if a high order bandpass filter with an irregular coupling topology is concerned.

The continuation method first constructs a set of polynomial system with coefficients that are functions of nonzero entries of a canonical coupling matrix, which are considered as parameters. By tracking the change of the unknowns, which are the entries of the orthogonal transformation matrix for reconfiguring the canonical coupling matrix to the one in the wanted coupling topology, with respect to changes in the parameters from the start system, the solution to the current system can be traced efficiently. The proposed algorithm adopts the general concept of the coefficient-parameter homotopy, a mathematic method for finding a specified solution from a known solution (the start system) in a different state. To overcome the difficulty of finding the start system with a solution in a different state of a filter with irregular coupling topology, a new concept called preconditioner is proposed for the first time. The preconditioner is physics-based and can greatly facilitate the continuation process for swiftly finding the coupling matrix that matches the physical realization.

The proposed theory will be validated through three practical examples, including one synthesis example and two model extraction examples for RAT. The first synthesis example is provided for illustrating the basic steps in applying the continuation approach with the proposed preconditioner. The rest of the two tuning examples concern RAT of 10-pole coaxial resonator filters with 6 or 7 TZs for 5G systems, demonstrating the effectiveness and efficiency of the proposed method for advanced industrial applications. The convergence of the solution and the bootstrapping approach for dealing with multiple irregular couplings are discussed in the examples.

II. PROPOSED RECONFIGURATION METHOD

The proposed reconfiguration method starts from establishing a set of polynomial equations according to the coupling topology concerned and the orthogonal transformation matrix to be sought. Unlike the existing numerical methods, whose starting point lacks legitimacy and the convergence is not guaranteed, the proposed method is based on a coefficient-parameter homotopy with a legitimate preconditioner for the start system, leading to fast convergence. The most attractive feature of the proposed method is its inheritability of a solution, which allows the solution at one tuning stage to serve as the start system of the subsequent tuning stage, being perfectly suitable for RAT of bandpass filters.

A. Parameterized Polynomial System

Following the formulation in [START_REF] Cameron | An exhaustive approach to the coupling matrix synthesis problem application to the design of high degree asymmetric filters[END_REF], n nonzero entries of the orthogonal transformation matrix Q are set to be the unknown vector = [ , ⋯ , ] . The multi-variable simultaneous polynomial equation system F(x, c(m)) = 0 is established as

    , , 0 , T T C i j i j T      Q Q I 0 M Q Q (1)
where M C is a coupling matrix in a canonical form, e.g. the folded form, T is the set of indices of coupling coefficients that are supposed to be zero in the target coupling matrix, and I is an identity matrix. The coefficient vector c is composed by the nonzero coupling elements m in M C , which is considered as the physical parameter vector of the current state. Obviously, the physical parameter vector m is different for different M C obtained from different filter responses, which results in different coefficients c for different polynomial systems. In this regard, the polynomial equation system F(x, c(m)) = 0 can be used to define the reconfiguration problem for different tuning states of the same filter.

B. Coefficient-Parameter Homotopy Continuation

The homotopy continuation is a primary computational method used in numerical algebraic geometry, in which a homotopy is formed between any two systems, and the isolated solutions of one are continued to those of the others. According to the polynomial system in (1), the coefficients of the polynomials defined by c(m) represent a system. Since c is a continuous function of m, a continuous path through the parameter space leads to a continuous evolution of the coefficients and multiple continuous paths for different solutions. A coefficient-parameter homotopy function can be defined as [START_REF] Sommese | Coefficient-parameter polynomial continuation[END_REF] 

     0 1 , , 1 H t F t t        x x c m m (2)
where m 1 is the parameter vector corresponding to a known solution with t = 1, which is called the start system, and m 0 is the parameter vector corresponding to the target system with t = 0. By decreasing t from 1 to 0, the homotopy function becomes the polynomial system F(x, c(m 0 )), in which x represents the vector containing the target unknowns to find. Normally, the solution x 1 to the start system is supposed to be known beforehand by solving the simultaneous polynomial equations. Usually, such a solution is obtained by an expensive numerical method.

The Euler-Newton path tracking algorithm [START_REF] Sommese | The Numerical Solution of Systems of Polynomials Arising in Engineering and Science[END_REF] illustrated in Fig. 1 is applied to track the solution to the polynomial system in (1) from t = 1 to 0. Taylor series expansion of the homotopy is applied by

        , ,
Higher-Order T , ,

erms t H t t H t H t t H t           x x x x x x x (3)
where

    , , H H t F       x x x c x x (4)
is the Jacobian matrix of H with respect to x and

        1 0 , , , t H H t t F t F         x c m x c m x ( 5
)
Suppose that the solution x i to the previous system at t i with respect to c i is known, or H(x i , t i ) = F(x i , c i ) = 0, an approximate solution i  x at t i+1 = t i + Δt can be predicted by setting H(x i +Δx,

t i +Δt) = 0 so that     1 , , i i t i i t t H H t      x x x x (6)
and the Euler prediction

i i     x x x (7)
Since ( , ) is not sufficiently small, the solution can be further updated by setting Δt = 0 in (3) at t i+1 . Then, the prediction ( 7) can be updated by

    1 1 1 , , i i i i t H t H         x x x x (8)
and the Newton correction

1 i i      x x x (9)
The pseudocode of the complete continuation process is given in Table I. By gradually changing t from 1 to 0, the solution corresponding to each t is found using the Euler-Newton method, and the expected solution corresponding to t = 0 can be eventually sought.

C. The Preconditioner for Start System

Suppose the coupling matrix M 1 of a response of the filter in the desired topology is known, the start system can be established by transforming M 1 to a chosen canonical form M C1 , which can always be done deterministically with transformation matrix Q defined by [START_REF] Cameron | General coupling matrix synthesis methods for Chebyshev filtering functions[END_REF] 1

T C  QM Q M (10) or 1 1 T C  Q M Q M with T  Q Q I (11)
in which the nonzero entries in Q compose solution x 1 and the nonzero entries in M C1 compose m 1 such that the start system of F(x 1 , c(m 1 )) = 0 is established. However, for most of the practical cases with irregular coupling topologies, such an M 1 is not easy to come by. Although the methods in [START_REF] Cameron | An exhaustive approach to the coupling matrix synthesis problem application to the design of high degree asymmetric filters[END_REF]- [START_REF] Pflüger | Coupling matrix extraction and reconfiguration using a generalized isospectral flow method[END_REF] can be used to find multiple possible solutions theoretically, the correct solution that matches the filter physical realization can only be picked out by exhaustively examining every solution, which is not practical 2) A set of parameters m 0 for the system to solve for RAT. In this work, a physics-based preconditioner for the start system is proposed.

It is very common in practice that an irregular coupling topology composed of a part with regular coupling topology and some additional weak irregular cross couplings. Physically, the filter response is dominated by the regular part whereas the intentionally introduced irregular cross couplings serve as "the icing on the cake". Therefore, the coupling matrix in the regular topology that aligns with the physical solution can be first reconfigured deterministically to construct a preconditioner of the start system. It will be shown through numerical examples that such a preconditioner can very well serve as the start system which effectively leads to the correct solution through the homotopy continuation process. When unintentional irregular cross couplings have to be introduced due to stray couplings, the solution corresponding to the part with regular coupling topology can be obtained by ignoring the far away TZs, leading to a physical preconditioner.

In dealing with an irregular coupling topology with multiple irregular cross couplings, in order to avoid multiple solution problem, a bootstrapping approach can be taken, by which the homotopy continuation process is applied multiple times with a progressive increase of irregular cross couplings. The details will be illustrated in the third example in the next section.

Note that this method is based on the hypothesis that the irregular coupling topology has a real solution for any lossless responses that comply with the shortest path rule [START_REF] Amari | On the maximum number of finite transmission zeros of coupled resonator filters with a given topology[END_REF]. If the irregular cross coupling is poorly located, there is a conceptual possibility that the topology has no physical solution. Nevertheless, an appropriate topology can eventually be obtained by trial-and-error since the failure of convergence of the proposed method indicates the bad choice of topology.

III. EXAMPLES AND CONVERGENCE OF SOLUTION

The proposed method is a general scheme not only for circuit model extraction problems but also for direct synthesis problems of bandpass filters with irregular topologies. The first illustrative example will be a synthesis problem of a 5-pole filter with an intentional irregular cross coupling, through which every detail of establishing a homotopy and the continuation process will be given. The second example is a practical 10-pole filter in a regular topology yet with an extra parasitic cross coupling, with which the accuracy of the method is demonstrated. The last example demonstrates how the proposed method is used in a RAT process of a 10-pole filter with more than one irregular couplings, in which the bootstrapping approach is applied.

A. Synthesis of a Five-Pole Filter with Three TZs

The synthesis of a five-pole filter that generates 3 TZs in the topology depicted in Fig. 2 is considered. Two of the three TZs are placed at -1.6j and 2.4j in the normalized lowpass frequency domain, which can be easily realized by two CTs. It's convenient to add a cross coupling between resonators 1 and 5 to create the third TZ at -8j to enhance the rejection in the far lower rejection band. For such a coupling topology there is no known rotation recipe for reconfiguration.

With the given TZs and specified 20 dB return loss level, the coupling matrix in the folded canonical form can be easily obtained and is denoted as M C . The first step to reconfigure M C to the coupling matrix, say M 0, in the desired topology is to establish a polynomial equation system. Basically, the orthogonal transformation matrix Q can be set as a blockwise diagonal matrix diag (1, Qx, 1) with Qx being an N × N block matrix whose elements are unknown [START_REF] Cameron | An exhaustive approach to the coupling matrix synthesis problem application to the design of high degree asymmetric filters[END_REF]. However, the set of the unknowns can usually be reduced to a less-redundant one by reducing the dimension of Qx such that Q = diag (I, Qx, I) in a trial-and-error way with symbolic calculation. In this example, the orthogonal transformation matrix Q can be set in the form of 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1

x x x x x x x x x                  Q (12)
where the unknown vector is x = [x 1 , x 2 , …, x 9 ]. Comparing Q T M C Q against the desired coupling topology, the homotopy function H is formulated according to (2) as

0 0 0 0 0 0 1 0 1 0 1 0 a x a x x a x x a x x a x x a x x a x x a x x a x x x x x x x x x a x x a x a x x x x x x x H x x x x x x x x x x x x x x                                                  (13) with a = t m 1 + (1 -t) m 0 and m = [M s,1 , M 1,1 , M 1,2 , M 2,2 , M 2,3 , M 3,3 , M 3,4 , M 4,4 , M 4,5 , M 5,5 , M 5,L , M 2,4 , M 2,5 , M 1,5 ]
. Subsequently, the Jacobian matrices with respect to x and t can be found according to (4) and [START_REF] Tamiazzo | An analytical technique for the synthesis of cascaded N-tuplets cross-coupled resonators microwave filters using matrix rotations[END_REF].

Assume the TZs 2.4j and -1.6j are assigned to CT (1-2-3) and CT (3-4-5), respectively. Then a physics-based approximation for M 1 is chosen to be the coupling matrix in the coupling topology without cross coupling M 1,5 , which can be synthesized as Fig. 3. Response of the start system M 1 (solid lines) and the target system M 0 (dashed lines) of the 5-pole filter.

                         M
By transforming M 1 to the folded form by a deterministic process as

C                          M
one can find the preconditioner Q for the start system, or   1 0.9594 0.2821 0 0.2530 0.8606 0.4421 0.1247 0.4241 0.8907

T    x
with parameter m 1 composed by coupling coefficients in M C1 .

Note that this is a degenerate system since it satisfies the polynomial equation system without cross coupling M 1,5 . The filter response of M 1 is shown in Fig. 3

C                          M
and consequently, the parameter m 0 , which represents the target system to solve. With the homotopy function and the preconditioner x 1 for the start system, the solution for the target system m 0 can be found by the coefficient-parameter homotopy continuation method as

                      

M

whose filter response is shown in Fig. 3 in dashed lines. As expected, higher rejection in the lower rejection band is obtained. By perturbing M 2,2 and M 4,4 in M 0 separately, it is confirmed that TZs corresponding to CT (1-2-3) and CT (3-4-5) are the same as those initially assigned, validating that M 0 is the correct solution.

It is interesting to look at the preconditioner from multiple solution point of view. By using a tool like Dedale-HF [START_REF] Seyfert | Dedale-HF[END_REF] one can compute all solutions associated to the irregular topology in Fig. 2. It is found that the total number (not counting the classical sign changes), which is the reduced order of the topology [START_REF] Cameron | An exhaustive approach to the coupling matrix synthesis problem application to the design of high degree asymmetric filters[END_REF], of solutions in this example is 2. When the third transmission zero at -8j is progressively sent to infinity, the two solutions converge to the two solutions corresponding to the two possible assignments of TZs to the two triplets. This explanation justifies the correspondence between the physics-based preconditioner in a regular coupling topology and the final solution in an irregular coupling topology. In other words, provided that the third TZ is sufficiently far away from the passband, the knowledge of the regular part of the coupling matrix can be used as a preconditioner to uniquely determine the final coupling matrix via the continuation process.

B. A 10-Pole Filter in Regular Topology with a Parasitic TZ

This example intends to demonstrate the effectiveness of the method in extracting the coupling matrix from the measured filter response. A coaxial combline filter in a multi-filter module that consists of 16 similar filters for a 5G base station product is shown in Fig. 4, the coupling topology of which shown in Fig. 5 contains one CT and two CQs in its physical realization. A parasitic coupling that is supposedly caused by the stray coupling between resonator 7 and the load, or M 7, L , is also presented in Fig. 5 in a dashed line. This parasitic coupling is chosen based on the shortest path rule [START_REF] Amari | On the maximum number of finite transmission zeros of coupled resonator filters with a given topology[END_REF] to accommodate the parasitic TZ mathematically. The center frequency of the filter f 0 = 2.591 GHz and the bandwidth BW = 194 MHz.

In this example, the polynomials for the filtering function are identified from the measured S-parameters using MVF with an extra TZ. With the polynomials, M C0 in the folded form can be easily obtained and is listed in Table II, with which the required polynomial homotopy system can be built according to [START_REF] Cameron | General coupling matrix synthesis methods for Chebyshev filtering functions[END_REF].

It can be easily found by perturbation that the desired TZ at -1.865 j is realized by CT (4-5-6); TZ pair at (1.136j, +1.237j) are controlled by CQ (1-2-3-4); and TZ pair at (-1.108j, -1.261j) are controlled by CQ (7-8-9-10). The coupling matrix M 1 in coupling topology without M 7, L can be easily synthesized with 20 dB return loss in the specified passband. The solution x 1 for the start system can be obtained by transforming M 1 to M C1 , which is in canonical form according to [START_REF] Atia | Synthesis of general topology multiple coupled resonator filters by optimization[END_REF]. Note that such a start system is a degenerate system since M 1 is not exactly in the desired coupling topology. However, the degenerate system can be found efficiently and can serve as a very good preconditioner.

The coupling matrix M C0 in canonical form is used as the parameters for the target system in the desired topology. By following the procedure given in Table I and using the preconditioner as the start system, the target system M 0 in the coupling topology of Fig. 5 can be obtained as listed in Table III, in which the imaginary part of the self-coupling terms represents the loss. Due to numerical errors in MVF caused by measurement noise and non-physical locations of spurious couplings to a certain extent, a very small amount of imaginary parts may occur in off-diagonal terms. Such small imaginary parts can be omitted without visibly affecting the response. As presented in Fig. 6, the filter response obtained by the reconfigured M 0 with the imaginary parts of inter-resonator couplings neglected matches the measured response very well.

In fact, a direct annihilating rotation sequence can be found for this special irregular coupling topology. After the conventional reconfiguration process for obtaining CQ (1-2-3-4) and CT (4-5-6), the final coupling matrix can be reconfigured by 3 more sequentially rotations, as listed in Table IV, for annihilating M 8,L , M 9,L , and M 8,10 , resulting in exactly the same coupling matrix as M 0 obtained by the homotopy continuation method. This fact further verifies the proposed preconditioner and the numerical continuation process.

Actually, the existence of coupling M 7,L changes the CQ closest to the load port to a quintet in arrow form, which is also a canonical block with a unique solution. Therefore, the degenerate system M 1 will evolve to the unique solution if the assignment of TZs to each cascaded block in M 1 is the same as that of M 0 . In fact, there is a one-to-one mapping in solutions between the irregular topology and any other topologies with which the missing cross coupling down the minimum path causes one TZ to vanish consequently. For RAT applications, this feature means that if a regular part of the coupling topology can be clearly identified the coupling matrix for the irregular topology will converge to the unique solution with the preconditioner.

The parasitic coupling M 7,L is artificially introduced to the circuit to accommodate the unintentional TZ. It's observed that the coupling value is relatively small and stable during the tuning process. A new golden template matrix is subsequently generated by incorporating the parasitic TZ. Then, the tuning process can be properly guided with the updated target coupling matrix.

C. A 10-Pole Filter with both Intentional and Unintentional Irregular Couplings

Another 10-pole filter is studied in this example, whose physical layout is shown in Fig. 7. This example will illustrate that the coefficient-parameter homotopy continuation can be efficiently applied to model extraction of a high order filter which is with both intentional and unintentional irregular couplings.

As illustrated in Fig. 8, the coupling topology consists of an intentional irregular cross coupling M 2,9 on top of the regular one CT and two CQs topology, as well as a parasitic coupling M 2,10 , which needs to be introduced to accommodate a parasitic TZ. Coupling M 2,9 is introduced conveniently to create an extra TZ to further suppress the rejection since resonators 2 and 9 are physically adjacent.

The extraction process will take two steps in a bootstrapping manner:

Step 1: a golden template coupling matrix that corresponds to the physical layout without considering the parasitic coupling M 2,10 , is synthesized using the proposed preconditioner and the continuation method. Assigning a TZ at -1.816j to CT(1-2-3), a pair of TZs (1.128j, 1.253j) to CQ(3-4-5-6), and a pair of TZs (-1.134j, -1.385j) to CQ(6-7-8-9), the preconditioner M 1 for synthesizing the template can be obtained based on the regular part of the filter layout without M 2,9 , and M 2,10 of course. The coupling matrix M C0 in folded form for the system with the additional TZ at about 3.5j due to cross coupling M 2,9 can be easily found by the traditional synthesis method. The synthesized response is equal-ripple in the passband. Similar to example A in this Section, the template coupling matrix M 0 corresponding to the physical layout can be obtained using the coefficient-parameter homotopy continuation, which is listed in Table V. The responses of the golden template matrix M 0 (in dashed lines) is compared to that of the preconditioner of M 1 (in solid lines) in Fig. 9. It is observed that with M 2,9 introduced, the suppression in the higher rejection band is increased.

Step 2: To accommodate the parasitic coupling M 2,10 , which is identified from the parasitic TZ in the measured response of the real filter, the golden template matrix M 0 is used as the preconditioner for reconfiguring the final coupling matrix from the measured data. Obviously, the final coupling matrix involves M 2,9 and M 2,10 . Since the additional TZ introduced by M 2,10 on top of those of the preconditioner is far away from the passband, applying the homotopy continuation process will lead to a unique solution. The bootstrapping approach can be applied to other reconfiguration problems, in which more than one irregular couplings exist. Step 2 can be applied to different tuning stages of a RAT process. Responses of reconfigured coupling matrices and the measured data of four typical tuning stages of the filter are superimposed in Figs. 10(a) -10(d), showing excellent pertinence of the extracted coupling matrix in each stage. Note that the golden template matrix needs to be constantly updated by taking into account the parasitic coupling in the tuning process.

It is this proposed reconfiguration method that enables a progressively automatic tuning process in a deterministic fashion. The reconfigured coupling matrices in the four tuning stages are listed in Table VI. The elapsed time for reconfiguration of each tuning stage using the proposed method with tens of iterations is about 0.75 seconds on average by a MATLAB [START_REF] Matlab | The Mathworks[END_REF] code on a PC computer.

IV. CONCLUSION

In this paper, an effective and efficient numerical approach for coupling matrix reconfiguration is proposed by introducing a legitimate and physics-based preconditioner in the parameter-coefficient continuation process. The preconditioner can greatly facilitate the continuation process to most of the practical filter problems, which involve irregular coupling topologies. By introducing the preconditioner, which can be easily obtained using classical filter synthesis theory, for the start system, the continuation process will lead to the unique solution corresponding to the given physical realization. The proposed approach is particularly useful not only for robot automatic tuning (RAT) but also for synthesizing advanced bandpass filters with irregular couplings, which are introduced either intentionally or unintentionally. The effectiveness and efficiency of the proposed approach have been validated and demonstrated through three practical examples, including a synthesis example, an extraction example, and a RAT of a 10-pole filter with 7 TZs. It is amazingly seen that the preconditioner can serve as a start system of the continuation

Fig. 1 .

 1 Fig. 1. Path tracking using prediction and correction steps. TABLE I PSEUDO-CODE OF THE COEFFICIENT-PARAMETER HOMOTOPY Input: The polynomial equation system F(x, c(m)) = 0 1) The start system, which is a set of known parameters m 1 and the solution x 1 for m 1 .

Fig. 4 .

 4 Fig. 4. A clipped photo of the coaxial combline filter (for example B) in a mult-filter module for 5G applications.

Fig. 5 .

 5 Fig. 5. Topology of the ten-pole filter in Fig. 4 with one CT, two CQs, and a parasitic coupling represented in the dashed line.

  in solid lines. The coupling matrix M C0 for the system with cross coupling M 1,5 in folded form can be easily found by the traditional synthesis method

Fig. 6 .

 6 Fig. 6. Measured response (solid line) and the recovered response (dashed line) corresponding to M 0 without the imaginary parts of inter-resonator couplings.

Fig. 7 .

 7 Fig. 7. A clipped photo of the coaxial combline filter (for example C) in a multi-filter module for 5G applications.

Fig. 8 .

 8 Fig. 8. Topology of a ten-pole filter in Fig. 7 with an irregular cross coupling and a parasitic coupling represented in the dashed line.

Fig. 9 .

 9 Fig. 9. Response of the preconditioner (solid lines) corresponding to M 1 and the synthesized golden template (dashed lines) corresponding to M 0 of example C.

Fig. 10 .

 10 Measured data in solid lines and responses by MVF in dashed lines of four tuning stages.

  Student Member, IEEE, Ke-Li Wu, Fellow, IEEE, and Fabien Seyfert

	A Direct Preconditioner for Coupling Matrix
	Reconfiguration of Bandpass Filters With
	Irregular Couplings Using Continuation Method
	Yan Zhang,

TABLE II COUPLING

 II MATRIX M C0 OF EXAMPLE B

	M S1	M 12	M 23	M 34	M 45	M 56
	1.0406 -0.0029j	0.7662 -0.0005j	0.6148 +0.0002j	0.5440 -0.0004j	0.3990	0.5681 -0.0004j
	M 67	M 78	M 89	M 9,10	M 10,L	M 57
	0.1916	0.5050	0.7218	0.8559	1.0426	0.0566
	+0.0005j	-0.0008j	+0.0001j	-0.0003j	-0.0018j	+0.0008j
	M 47	M 48	M 38	M 39	M 29
	-0.5232	0.1237 +0.0008j	0.0844	-0.0312 -0.0003j	-0.0009 +0.0001j
	M SS	M 11	M 22	M 33	M 44	M 55
	-0.0016	0.4309	-0.5044	-0.2215	-0.4212	-0.1097
	-0.0142j	-0.0134j	-0.0073j	-0.0069j	-0.0075j	-0.0073j
	M 66	M 77	M 88	M 99	M 10,10	M LL
	0.6160	0.1665	-0.0512	0.0358	0.0152	0.0043
	-0.0082j	-0.0061j	-0.0068j	-0.0074j	-0.0196j	-0.0218j

TABLE III COUPLING

 III MATRIX M 0 OF EXAMPLE B

	M S1	M 12	M 14	M 23	M 24	M 34
	1.0406 -0.0029j	0.6601 -0.0009j	0.3889 +0.0005j	0.1533 +0.0002j	0.4080 -0.0007j	0.2240
	M 45	M 46	M 56	M 67	M 78	M 79
	0.5381	-0.1633	0.4533	0.5410	0.2607	-0.4610
	+0.0002j	+0.0005j	-0.0001j	-0.0006j	-0.0015j	+0.0024j
	M 7,10	M 7L	M 89	M 9,10	M 10,L
	0.3834	0.0170	0.1785	0.7751	1.0425
	+0.0027j	-0.0011j	-0.0016j	-0.0024j	-0.0018j
	M SS	M 11	M 22	M 33	M 44	M 55
	-0.0016	0.4309	-1.0311	-0.3959	-0.3722	0.1547
	-0.0142j	-0.0134j	-0.0082j	-0.0069j	-0.0069j	-0.0083j
	M 66	M 77	M 88	M 99	M 10,10	M LL
	-0.0757	-0.2720	0.9430	0.5717	0.0028	0.0043
	-0.0081j	-0.0104j	-0.0066j	-0.0028j	-0.0189j	-0.0218j

TABLE IV ROTATION

 IV SEQUENCE TO RECONFIGURE THE COUPLING MATRIX Rotation sequence

		Elements to be	Pivot
		annihilated	[i, j]
	1	[8, L]	[8, 9]
	2	[9, L]	[9, 10]
	3	[8, 10]	[8, 9]

TABLE V

 V

				COUPLING MATRIX M 0 FOR THE GOLDEN TEMPLATE MATRIX OF EXAMPLE C
		S	1	2	3	4	5	6	7	8	9	10	L
	S	0	0.9809	0	0	0	0	0	0	0	0	0	0
	1	0.9809	0	0.7999	-0.0937	0	0	0	0	0	0	0	0
	2	0	0.7999	0.1344	0.5723	0	0	0	0	0	-0.0034	0	0
	3	0	-0.0937	0.5723	-0.0179	0.2317	0.4392	0.2075	0	0	0	0	0
	4	0	0	0	0.2317	-0.9201	0.1449	0	0	0	0	0	0
	5	0	0	0	0.4392	0.1449	-0.3598	0.4843	0	0	0	0	0
	6	0	0	0	0.2075	0	0.4843	-0.018	0.2376	-0.4298	0.217	0	0
	7	0	0	0	0	0	0	0.2376	0.9064	0.1669	0	0	0
	8	0	0	0	0	0	0	-0.4298	0.1669	0.3512	0.5355	0	0
	9	0	0	-0.0034	0	0	0	0.217	0	0.5355	0	0.8053	0
	10	0	0	0	0	0	0	0	0	0	0.8053	0	0.9809
	L	0	0	0	0	0	0	0	0	0	0	0.9809	0

TABLE VI REAL

 VI PARTS OF COUPLING MATRICES OF THE FOUR TUNING STAGES IN FIGS. 10 (A)-(D) very robustly. It is believed that the approach can play an important role in future RAT of advanced microwave bandpass filters.

		Stage 1	Stage 2	Stage 3	Stage 4
	M SS	-0.0009	-0.0014	-0.0012	-0.0014
	M 11	-0.0259	-0.0219	0.0017	0.0171
	M 22	-0.0341	0.2967	0.2768	0.3683
	M 33	0.1567	0.1242	-0.0526	-0.0420
	M 44	-0.6493	-0.8493	-0.9499	-0.9607
	M 55	0.1318	-0.2775	-0.4011	-0.4142
	M 66	-0.2617	0.0865	-0.0210	-0.0134
	M 77	0.5715	0.8695	0.9461	0.9337
	M 88	0.3370	0.6802	0.4027	0.3700
	M 99	0.0136	0.0120	0.0135	0.0108
	M 10,10	0.0492	0.0547	0.0044	0.0204
	M LL	0.0025	0.0071	-0.0074	0.0132
	M S1	1.0594	1.0592	1.0589	1.0586
	M 12	0.7882	0.8326	0.8431	0.8430
	M 23	0.5091	0.5034	0.5410	0.5740
	M 34	0.2235	0.2187	0.2149	0.2165
	M 45	0.1680	0.1654	0.1652	0.1265
	M 56	0.3889	0.4249	0.4253	0.4854
	M 67	0.1913	0.1833	0.1953	0.2387
	M 78	0.1760	0.1802	0.1693	0.1677
	M 89	0.5581	0.5708	0.5538	0.5537
	M 9,10	0.8971	0.8967	0.8980	0.8564
	M 10,L	1.0069	1.0068	1.0047	1.0167
	M 13	-0.2369	-0.2362	-0.2043	-0.2209
	M 35	0.4497	0.4525	0.4559	0.4523
	M 36	0.1859	0.1958	0.1833	0.2304
	M 68	-0.4780	-0.4787	-0.4708	-0.4533
	M 69	0.2580	0.2431	0.2717	0.2330
	M 29	-0.0116	-0.0107	-0.0027	-0.0035
	M 2,10	-0.0007	-0.0007	0.0009	0.0009

process


J J Verbeek

N Vlassis

B Kröse

A k-segments algorithm for finding principal curves

Keywords: dimension reduction, feature extraction, polygonal line, principal curve, unsupervised learning

We propose an incremental method to find principal curves. Line segments are fitted and connected to form polygonal lines. New segments are inserted until a performance criterion is met. Experimental results illustrate the performance of the method compared to other existing approaches.

Introduction

Principal curves are the nonlinear generalization of principal components. They give a summarization of the data in terms of a 1-d space non-linearly embedded in the data space. Intuitively, a principal curve 'passes through the middle of the (curved) data cloud'. Applications include dimension reduction for feature extraction, visualization and lossy data compression. What concerns the first application, the curve can be used to obtain an ordering of the data by projecting the data to the curve. An example of such a use of principal curves is the ordering of ecological species abundance data [Deáth, 1999]. Applications of the proposed principal curve algortihm can also be envisaged in fitting and recognizing complex shapes in images. Several definitions of principal curves have been proposed in the literature. One of the earliest definitions is based on 'self-consistency' [START_REF] Hastie | [END_REF], i.e. the curve should coincide at each position with the expected value of the data projecting to that position. Another approach [Kégl et al., 2000] is to define principal curves of length l as curves of length l that achieve the minimum expected squared distance from points to their projection on the curve. The Polygonal Line Algorithm (PLA) finds principal curves in terms of the latter definition. A more probabilistic approach [Tibshirani, 1992] defines principal curves as curves minimizing a penalized log-likelihood measure. The data is modeled by a mixture of the form:

p(x) = l 0 p(x|t)p(t)dt ( 1
)
where t is a latent variable distributed on an arc-length parameterized curve of length l. p(x|t) is, for example, a spherical Gaussian modeling the noise located on point t of the curve.

In addition to the above mentioned principal curve algorithms, several other methods can be used to tackle this and closely related problems. The Generative Topographic Mapping (GTM) [Bishop et al., 1998], finds a nonlinear curve embedded in the data space. Along this curve a mixture of Gaussian kernels is placed to generate a distribution on the data space. Also some vector quantization techniques can be used to find approximations to principal curves, e.g. Self

Organizing Maps (SOM) [Kohonen, 1995] and Growing Cell Structures (GCS) [Fritzke, 1994].

These methods incorporate a topology among the prototype-vectors. The original space is then mapped to a discrete 1-d space instead of a continuous 1-d space.

Several problems exist with the afore mentioned methods. A common feature of these is that they consist of a combination of 'local models' that are related by a fixed topology. When the data is concentrated around a highly curved or self-intersecting curve, these methods exhibit poor performance. This is due to the fixed topology among the local models and to bad initialization. Also, often one does not know a-priori how many 'local models' one needs and one has to make an (educated) guess for the number of local models. On the other hand, the definition of [START_REF] Hastie | [END_REF]] explicitly excludes self-intersecting curves as principal curves. In general PLA and Tibshirani's method also fail due to bad initialization in cases where the data is concentrated around highly curved or self-intersecting curves. See Figure 1 for an illustration (for PLA, the dashed curve passing through the data is the generating curve).

In this paper, we use a probabilistic setting to find principal curves by means of maximizing log-likelihood, resembling the method of [Tibshirani, 1992]. We assume the data is corrupted by some noise, therefore we do not require the curve to fit the data exactly. We use the model (1), where p(t) is uniform along the curve and p(x|t) is a spherical Gaussian located on point t of the curve, with constant variance over all t. The variance is a smoothing parameter, to be set by the user. We restrict the considered curves to polygonal lines (PLs). Our method avoids many of the poor solutions mentioned above (and shown in Figure 1) by fitting local linear models without any topological constraints. To enable an efficient search among PLs with a varying number of segments and to avoid local minima, we use an incremental method that considers PLs consisting of an increasing number of segments. Hence, it starts with models that are relatively easy to fit and gradually increases their complexity. Overview of the method. By means of local Principal Component Analysis (PCA) we find local linear models. The local linear models are fitted by a (two-step) iterative update scheme.

After fitting, the local linear models are combined to form a PL by means of a heuristic search.

Next, if appropriate, a new local linear model is added. The location of the new local linear model is found by means of a global search over the data set. Among all PLs encountered by the algorithm, the PL that maximizes an approximation of the log-likelihood of the data is selected.

Overview of the paper. In the next sections we discuss the algorithm in detail. Section 2 is concerned with how we fit the local linear models. In Section 3 we describe the method to find 'good' PLs given a set of line segments. In Section 4 the objective function is considered. We end the paper with a discussion.

A k-segments algorithm

In order to find a 'good' PL for the data, we first search for a set of k 'good' line segments. In the first subsection we discuss how we extend the k-means algorithm [Gersho and Gray, 1992] to an algorithm for finding k lines. In Section 2.2 we adapt the k-lines algorithm to find line segments that can be used for PL construction. In the last subsection we describe a method to insert a new segment given k segments and a data set.

Extending k-means to k-lines.

The Generalized Lloyd Algorithm or 'k-means' is a well-known vector quantization method. To extend this method to k-lines, we start with some definitions. A line s is defined as: s = {s(t)|t ∈ IR}, where s(t) = c + ut. The distance from a point x to a line s is defined as:

d(x, s) = inf t∈IR s(t) -x .
(2)

Let X n be a set of n samples from IR d . Define the Voronoi Regions (VRs) [Gersho and Gray, 1992] V 1 , . . . , V k as:

V i = {x ∈ X n | i = arg min j d(x, s j )}. (3)
Hence, V i contains all data points for which the i th line is the closest, see Figure 2 for an illustration.

Analogue to k-means, the goal is to find k lines s 1 , . . . , s k that minimize the total squared distance of all points to their closest line:

k i=1 x∈Vi d(x, s i ) 2 . ( 4
)
To find lines that are local optima of (4) we modify the k-means algorithm slightly: Start with random orientations and locations of the k lines. Next, alternate between the following two steps until convergence:

1. determine the VRs.

replace the lines by the first Principal Component (PC) of their VR. 1

To prove that this algorithm converges, we make two observations:

• The objective function (4) decreases both in step 1 and step 2. For step 1 this follows from the definition of the VRs. For step 2 this follows from the reconstruction error minimization property of the first PC [Ripley, 1996].

• Due to the finite cardinality of X n , there are only a finite number of distinct (Voronoi) partitions of X n .

From these two observations it follows trivially that the algorithm achieves a local optimum for (4) in a finite number of steps.

2.2 From k-lines to k-segments.

Since we want to construct a polygonal line we actually do not search for lines but for line segments instead. We therefore replace step 2 of the algorithm. Instead of using the first PC, we

1 More precisely: we first translate the data in the VR to obtain zero mean and then take the line along the first PC of the translated data. use a segment of it. For example, we could use the shortest segment of the first PC such that all orthogonal projections to the first PC of the points in the VR are included in the segment. We adjust the definition of the VR straightforwardly by using the distances to the segments, i.e. t is restricted to a bounded interval in (2). Using these segments the value of the objective function is not altered. However, two observations suggested another approach:

V 1 V 0 V 2
• The objective function ( 4) is an easy to compute substitute for the log likelihood of the data given the segments, see Section 4. The length of the segments is not reflected in (4). As a result relatively long segments are found which result in sub-optimal PLs.

• In applications of the algorithm we sometimes observed poor performance. This was due to local optima of (4) which could be avoided by removing the restriction that the segments have to include all projections of the points in the VR to the first PC.

To include explicit optimization of the segment length would be computationally very costly.

Therefore, we settled for a heuristic that works well in practice. Experimentally, we found that in general good performance is obtained if we use segments of the first PC that are cut off at 3σ/2 from the centroid of the VR, where σ 2 is the variance along the first PC. To maintain the convergence property, we check whether ( 4) is actually decreased by the new segment. If (4) is not decreased, we use the segment that includes all projections to the first PC to obtain guaranteed decrease of (4).

Incremental k-segments.

In most applications we do not know in advance how many segments are 'optimal' to model the data. Therefore, it makes sense to try (many) different numbers of segments and use some criterion to select the best number. Since running the algorithm for many different numbers of segments is time consuming, we propose to use an incremental strategy. We start with k = 1 segment and use the search method described in the previous subsections to find optimal segments. Next, we insert a new segment and optimize the k + 1 segments again. The scheme is repeated until some criterion is met (i.e. a predefined maximum of lines is reached or some performance criterion is met).

To determine where to insert the new segment, we first compute for each x i in our data set the decrease of (4) if we would place a zero-length segment on x i . Below we explain why we use this criterion. With a zero-length segment we just mean a point c, hence the distance function

(2) reduces in this case to d(x, s) = xc . The VRs of the zero-length segments are defined as before, by (3), see also Figure 2. To avoid obtaining segments representing too few points, we only consider those zero-length segments for which the cardinality of the corresponding VR is at least three. Let V k+1 be the VR of the zero-length segment maximizing the decrease in (4). Then we insert a segment along the first PC of the data in V k+1 which is cut off at 3σ/2 on each side of its mean, where σ 2 is the variance of the data in V k+1 when projected orthogonally to the first principal component of the data in V k+1 . After the insertion of the new segment, again we are ready to do step-wise optimization of the segments. Note that the zero-length segments are not inserted themselves, but only used to select a good position to insert the new segment.

The decrease in (4) due to a insertion of a zero-length segment provides a lower bound on the decrease due to the segment insertion that is actually performed. To see this, let V k+1 denote the VR corresponding to the zero-length segment at x i . It is a well known fact [Gersho and Gray, 1992] that for a given finite set of points S ⊂ IR d the mean m of those points minimizes the squared distance function: m = arg min µ∈IR d x∈S x -µ 2 , hence:

x∈S x -m 2 ≤ x∈S x -x i 2 .
Since the new segment s includes m, the total square distance to the segment is at most the total square distance to m i.e. x∈S d(x, s) 2 ≤ x∈S xm 2 . Hence, we have the lower bound. Note that the search for the minimizing x i can be implemented efficiently as described in Appendix A.

Combining segments into a polygonal line

Until now we have introduced an incremental method to find a set of line segments. The outset, however, was to find a PL maximizing the log-likelihood of the data. In this section we discuss how to link the segments together to form a PL Since there are in general many ways to link the segments together, we settle for a greedy strategy. First we generate an initial PL and next we make step-wise improvements on the initial PL.

To achieve a fast algorithm to find the PL we do not consider the data in the construction of the PL. We define a fully connected graph G = (V, E), where the set of vertices V consists of the 2k end-points of the k segments. Also, define a set of edges A ⊂ E which contains all edges that V V' V V' correspond to the segments. A sequence of edges

{(v 0 , v 1 ), (v 1 , v 2 ), . . . , (v m-2 , v m-1 ), (v m-1 , v m )}
in which all edges are distinct is called a 'path' [Wilson, 1972].

A path is 'open' if v 0 = v m . An
open path that passes through every vertex in the graph exactly once is called a 'Hamiltonian path' (HP). Note that we can consider a HP as a set P ⊂ E. We want to find the HP P minimizing the total cost of the path, under the constraint:

A ⊂ P ⊂ E. 2
The cost of a path P is defined as l(P ) + λa(P ), with 0 ≤ λ ∈ IR is a parameter to be set by the user. The term l(P ) denotes the length of the path, defined as the sum of the lengths of the edges in P . The length of an edge e = (v i , v j ) is taken as the Euclidean distance between its vertices: l(e) = v i -v j . The second term, a(P ), is a penalty term equal to the sum of the angles between adjacent edges. The parameter λ controls the trade-off between preferring short paths and paths that do not contain sharp turns. We introduced the λa(P ) term to implement a preference for 'smooth' (not having sharp turns) curves. The smaller λ, the smaller the preference for smooth curves. Especially in cases where the generating curve is self-crossing (c.f. Figure 5) setting λ > 0 is important to find the right PL.

To construct an initial PL we use the greedy strategy outlined below. We call a HP on a subset of V a sub-HP. We start with the original k segments as k sub-HPs. At each step we connect two sub-HPs with an edge e. Note that the total cost of a (sub-) HP consisting of two sub-HPs P i and P j linked together by an edge e is the sum of the costs of each sub-HP plus l(e) plus an angle penalty a(e). Figure 3 illustrates the angle penalty a(e) = α + β incurred by an edge e = (v i , v j ) that connects two sub-HPs P i and P j . We assign to each edge e ∈ (E -A) cost c(e) = l(e) + λa(e).

The construction algorithm inserts at each step the edge that minimizes c(e) over all edges that connect two sub-paths. Summarizing, the procedure is as follows:

1. Start with k sub-HPs defined by A.

2. While there are at least two sub-HPs 3. Join those two sub-HPs P i and P j (i = j) by edge e ∈ (E -A) such that e minimizes c(e) over all edges connecting two distinct sub-HPs.

Due to the greedy nature of the construction discussed above, it might well happen that we can obtain a HP with lower cost by some simple modifications of the initial HP. To find such improvements we use a simple result from graph theory [START_REF] Lawler | [END_REF]: finding the optimal (in terms of cost) HP for k vertices, can be expressed as finding the optimal solution for a k+1 cities traveling salesman problem (TSP). We use the 2-opt TSP optimization scheme [START_REF] Lawler | [END_REF] to improve our initial HP.

Objective function

Using the incremental method described above, we obtain a sequence of PLs with an increasing number of segments. As stated in the introduction, we search for the PL that maximizes the log-likelihood of the data.

Consider a PL of length l as a continuous arc-length parameterized 1-d latent variable t embedded in IR d , where the embedding is given by f : [0, l] → IR d . For simplicity, we assume a uniform distribution p(t) = 1/l on t here for t ∈ [0, l]. Furthermore, let p(x|t) be a spherical Gaussian distribution with mean point t on the PL and covariance matrix σ 2 I, where I is the d × d identity matrix. It turns out that, for a latent variable distributing uniformly over a line segment s of length l σ, the negative log-likelihood for a point x can be roughly approximated by

log l + d(s, x) 2 /(2σ 2 ) + c, (5)
where c is some constant dependent on σ. This can be seen as follows:

-log p(x) = -log t∈[0,l] p(x | t)p(t)dt = log(l) + c 1 -log t∈[0,l] exp( x -f (t) 2 2σ 2 )dt. (6)
We may write xf (t) 2 as the sum

d ⊥ (f (t), x) 2 + d (f (t), x) 2 , see Figure 4. Since d ⊥ (f (t), x) = d ⊥ (s, x
) is constant for different t, we can write the last term in (6) as:

d ⊥ (s, x) 2 /(2σ 2 ) -log t∈[0,l] exp(-d (f (t), x) 2 )dt. (7) Let d (s, x) = inf t∈[0,l] {d (f (t),
x)}, then the last integral can be approximated roughly by

c 2 exp(-d (s, x) 2 /(2σ 2 )). Note that d(s, x) 2 = d (s, x) 2 + d ⊥ (s, x) 2 , this leads to equation (5).
We neglect the effect of higher density occurring at one side of the PL at places where different non-parallel segments are connected. Hence, we approximate the total log-likelihood of the data as:

n log l + k i=1 x∈Vi d(s i , x) 2 /(2σ 2 ), ( 8
)
where l is the total length of the PL.

In practice, if we take the objective function as a function of k, the number of segments, the global minimum often is the first local minimum of the objective function. So, a simple stopping

f (t) x d(s, x) d ⊥ (s, x) s d (f (t), x) d (s, x) d(f (t), x)
Figure 4: Illustration of the distances between a point x a segment s and a point f (t) on s.

-

1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -1.5 -1 -0.5 0 0.5 1 1.5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -1.5 -1 -0.5 0 0.5 1 1.5
Figure 5: Results on a synthetic data set for our method (left) and GTM (right). strategy could be to keep inserting segments until the objective function (8) reaches its first minimum. It is also possible to keep inserting segments until some limit k max on k is reached. A variety of other stopping criteria is possible.

Discussion

To test the algorithm presented above, we conducted experiments on several artificial data sets.

We used 2-d and 3-d data spaces, crossing curves, non-crossing curves along which we generated data, different amounts of noise and different numbers of training examples. The results are promising, see Figure 5 and 6 for an illustration (the thick segments in the plots are the fitted segments, the thin segments are inserted by the PL construction). Due to limited space we cannot present exhaustive experimental results but just some illustrations here. The good performance is due to (i) the flexibility of the method: the segments are fitted without being connected, and (ii) its incremental nature: segments are added one by one to avoid local optima and (iii) the 'optimal' number of segments is determined automatically. In Figure 6 we illustrate, on the same data as in Figure 1, several stages of the principal curve fitting.

A difference between our method and other principal curve algorithms is the absence of a curvature penalty in our objective function (8). Although a curvature penalty is used in the con-struction of the curves, the algorithm is not too sensitive to changes in its setting and it may even be omitted when modeling non-crossing curves. Regularization of the model complexity is instead introduced in ( 8) by (i) modeling the distribution of the latent variable by a uniform distribution and (ii) assuming a given level of noise all along the curve. This results in a regularization parameter σ with a clear physical interpretation. The essence is that we avoid (i) distributions on the latent variable that have only finite support 3 and (ii) arbitrarily small noise levels. Experimentally we found that the applicability of the method is not limited to cases where the naive assumption of a uniform distributed latent variable holds. Furthermore, it can be readily replaced by less rigid assumptions on the distribution of the latent variable, for example by a segment-wise uniform distribution on the latent variable.

Note that when one wants to fit other structures than curves, for example closed paths or branching structures, one only needs to replace the search procedure described in Section 3. In the case of closed paths the modifications needed are trivial since the same construction and optimization methods can be used.

Our algorithm has a running time O(kn 2 ) that is dominated by the insertion procedure, see Appendix A. For each of the k insertions the allocation takes O(n 2 ) operations and hence takes O(kn 2 ) in total. Constructing and optimizing all the PLs takes in total O(k 3 ) (taking the maximal number of iterations of 2-opt fixed). Computing the VRs in the optimization procedure after each segment allocation takes O(nk 2 ) operations in total. The total number of operations needed by the determination of the principal components is O(kn) (taking the dimensionality of the data space fixed). Since in general k n we arrive at a running time of O(kn 2 ). In applications involving large amounts of data the quadratic scaling with the amount of data may be problematic. In such cases more sophisticated and faster insertion methods may be used to allow for a O(k 2 n) run time. For example only a small amount of candidate insertion locations may be considered using a method analogue to the one used in [Verbeek et al., 2001].

Some extensions on the presented method suggest themselves:

Smooth curves In some situations there may be a preference for smooth curves over PLs. In such cases one can use the found PL to assign a latent variable value to each datum (implying an ordering on the data). Subsequently, one may use some regression method to find a 'good' smooth curve. The found PL may also serve as an initialization of one of the other methods mentioned that find smooth curves.

Obviate variance estimation The use of an objective function that does not need to be tuned by the user would increase usefulness of the method. We plan to use the Minimum Description

3 The distribution on the latent variable maximizing likelihood is a distribution with support only at the points on the curve where data projects [Tibshirani, 1992]. Soft version Currently we are investigating a 'soft' version of the presented method in which the 'hard' Voronoi partitions are replaced by a probabilistic weighting scheme. This way we hope to obtain better performance on very noisy data.

Software implementing the k-segments algorithm in Matlab can be obtained from the homepage of the first author: http://www.science.uva.nl/ ∼jverbeek.

A Efficient search for allocation of new segment

At the start of the algorithm we compute the pairwise squared distances between all points in the data set, denote the resulting symmetric n × n matrix by D. In the optimization of the segments, we determine their VRs at each step. To do so we compute distances between all points and all segments, we store for each point x i the squared distance d V R i to the closest segment. Let

D V R = (d V R 1 , d V R 2 , . . . , d V R n ) and VD = [D V R . . . D V R
] be a square n × n matrix. If A is a matrix and a ∈ IR, then let M = max(A, a) denote the matrix M where M i,j = max(A i,j , a). Let G = max(VD -D, 0). Then G i,j equals the decrease in squared distance for x i to its closest segment if we insert the zero-length segment at x j . Let 1 denote the 1 × n matrix [1 1 . . . 1]. Then, arg max{1G} gives us the index of the datum maximizing the lower bound.

Figure 1 :

 1 Figure 1: A spiral with some noise. Results for (left to right, top to bottom): GTM, PLA, GCS and SOM.

Figure 2 :

 2 Figure 2: An example where the data (open dots) is partitioned in three Voronoi Regions for two normal segments (V 1 and V 2 ) and one (V 0 ) 'zero-length' segment (black dot).

Figure 3 :

 3 Figure3: Connecting two sub-paths: the angle penalty of (v i , v j ) is the sum of the angles between an edge and the adjacent edges, i.e. α + β.

Figure 6 :

 6 Figure 6: Some of the PLs encountered (using 4, 6 and 12 segments respectively) by the algorithm when fitting data distributed along a spiral.

There always exists a HP for G that contains A since G is fully connected and A connects only pairs of vertices.

* This research is supported by the Technology Foundation STW, applied science division of NWO and the technology programme of the Ministry of Economic Affairs.


J J Verbeek

N Vlassis

B Kröse

Coordinating Principal Component Analyzers

Mixtures of Principal Component Analyzers can be used to model high dimensional data that lie on or near a low dimensional manifold. By linearly mapping the PCA subspaces to one global low dimensional space, we obtain a 'global' low dimensional coordinate system for the data. As shown by Roweis et al., ensuring consistent global low-dimensional coordinates for the data can be expressed as a penalized likelihood optimization problem. We show that a restricted form of the Mixtures of Probabilistic PCA model allows for a more efficient algorithm. Experimental results are provided to illustrate the viability method.

Introduction

With increasing sensor capabilities, powerful feature extraction methods are becoming increasingly important. Consider a robot sensing its environment with a camera yielding a stream of 100 × 100 pixel images. The observations made by the robot often have a much lower intrinsic dimension than the 10.000 dimensions the pixels provide. If we assume a fixed environment, and a robot that can rotate around its axis and translate through a room, the intrinsic dimension is only three. Linear feature extraction techniques are able to do a fair compression of the signal by mapping it to a much lower dimensional space. However, only in very few special cases the manifold on which the signal is generated is a linear subspace of the sensor space. This clearly limits the use of linear techniques and suggests to use non-linear feature extraction techniques.

Mixtures of Factor Analyzers (MFA) [START_REF] Ghahramani | The EM Algorithm for Mixtures of Factor Analyzers[END_REF] can be used to model such non-linear data manifolds. This model provides local linear mappings between local latent spaces and the data space. However, the local latent spaces are not compatible with each other, i.e. the coordinate systems of neighboring factor analyzers might be completely differently oriented. Hence, if we move through the data space from one factor analyzer to the next we cannot predict the latent coordinate of a data point on the one factor analyzer from the latent coordinate on the other factor analyzer.

Recently, a model was proposed that integrates the local linear models into a global latent space, allowing for mapping back and forth between the global latent and the data-space [START_REF] Roweis | Global coordination of local linear models[END_REF]. The idea is that there is a linear map for each factor analyzer between the data-space and the global latent space. The model, which is fitted by maximizing penalized log-likelihood with an algorithm closely related to the Expectation-Maximization (EM) algorithm [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF], is discussed in the next section. In Section 3, we show how we can reduce the number of parameters to be estimated and simplify the algorithm of [START_REF] Roweis | Global coordination of local linear models[END_REF]. These simplifications remove the iterative procedure from the M-step and remove all matrix inversions from the algorithm. The price we pay is that the covariance matrices of the Gaussian densities we use are more restricted, as discussed in the same section. Experimental results are given in Section 4. A discussion and conclusions are provided in Section 5.

The density model

To model the data density in the high dimensional space we use mixtures of a restricted type of Gaussian densities. The mixture is formed as a weighted sum of its component densities, which are indexed by s. The mixing weight and mean of each component are given by respectively p s and µ s . The covariance matrices of the Gaussian densities are constrained to be of the form:

C = σ 2 (I D + ρΛΛ ), Λ Λ = I d , ρ > 0 (1)
where D and d are respectively the dimension of the high-dimensional/data-space and the low-dimensional/latent space. We use I d to denote the d-dimensional identity matrix. The d columns of Λ, in factor analysis known as the loading matrix, are D-dimensional vectors spanning the local PCA. Directions within the PCA subspace have variance σ 2 (1 + ρ), other directions have σ 2 variance. This is as the Mixture of Probabilistic Principal Component Analyzers (MPPCA) model [START_REF] Tipping | Mixtures of probabilistic principal component analysers[END_REF], with the difference that here we do not only have isotropic noise outside the subspaces but also isotropic variance inside the subspaces. We use this density model to allow for convenient solutions later. In [START_REF] Verbeek | Procrustes Analysis to Coordinate Mixtures of Probabilistic Principal Component Analyzers[END_REF] we provide a Generalized EM algorithm to find maximum likelihood solutions for this model. The same model can be rephrased using hidden variables z, which we use to denote 'internal' coordinates of the subspaces. We scale the coordinates z such that: p(z | s) = N (z; 0, I d ). The internal coordinates allow us to clearly express the link to the global latent space, for which we denote coordinates with g. All mixture components have their own linear mapping to the global space, given by a translation κ and a matrix A, i.e. p(g | z, s) = δ(κ s + A s z), where δ(•) denotes the distribution with mass 1 at the argument. The generative model reads:

p(x | z, s) = N (x; µ s + √ ρ s σ s Λ s z, σ 2 s I D ), (2) p(x) = s p s N (x; µ s , σ 2 s (I D + ρ s Λ s Λ s )), p(g) = s p s N (g; κ s , A s A s ). ( 3
)
We put an extra constraint on the projection matrices:

A s = α s σ s √ ρ s R s , R s R s = I d , α s > 0, (4)
hence R s implements only rotations plus reflections.

Note that the model assumes that locally there is a linear correspondence between the data space and the latent space. It follows that the densities p(g | x, s) and p(x | g, s) are Gaussian and hence both p(g | x) and p(x | g) are mixtures of Gaussian densities. In the next section we discuss how this density model allows for an efficient learning scheme, as compared to the expressive but expensive MFA model proposed in [START_REF] Roweis | Global coordination of local linear models[END_REF].

The learning algorithm

The goal is, given observable data {x n }, to find a good density model in the data-space and mappings {A s , κ s } that give rise to 'consistent' estimates for the hidden {g n }. With consistent we mean that if a point x in the data-space is well modeled by two PCA's, then the corresponding estimates for its latent coordinate g should be close to each other, i.e. the subspaces should 'agree' on the corresponding g.

Objective Function: To measure the level of agreement, one can consider for all data points how uni-modal the distribution p(g | x) is. This idea was also used in [START_REF] Vlassis | Supervised dimension reduction of intrinsically low-dimensional data[END_REF]. There the goal was to find a global linear low-dimensional projection of supervised data, that preserves the manifold structure of the data. In [START_REF] Roweis | Global coordination of local linear models[END_REF] it is shown how the double objective of likelihood and uni-modality can be implemented as a penalized log-likelihood optimization problem. Let Q

(g | x n ) = N (g; g n , Σ n ) a Gaussian approximation of the mixture p(g | x n ) and Q(s | x n ) = q ns . We define: Q(g, s | x n ) = Q(s | x n )Q(g | x n ). (5)
As a measure of uni-modality we can use a sum of Kullback-Leibler divergences:

ns dg Q(g, s | x n ) log Q(g, s | x n ) p(g, s | x n ) = n D KL (q ns p ns ) + ns q ns D ns where D ns = D KL (Q(g | x n ) p(g | x n , s)
) and p ns = p(s | x n ). The total objective function, combining log-likelihood and the penalty term, then becomes:

Φ = n log p(x n ) -D KL ({q ns } {p ns }) - s q ns D ns (6) = ns dg Q(g, s | x n ) -log Q(g, s | x n ) + log p(x n , g, s) . (7)
The objective corresponds to a constrained EM procedure, c.f. [START_REF] Verbeek | The Generative Self-Organizing Map: A Probabilistic Generalization of Kohonen's SOM[END_REF] where the same idea is used to derive a probabilistic version of Kohonen's Self-Organizing Map [START_REF] Kohonen | Self-Organizing Maps[END_REF]. Our density model differs with that of [START_REF] Roweis | Global coordination of local linear models[END_REF] in two aspects: (i) we use an isotropic noise model outside the subspaces (as opposed to diagonal covariance matrix) and (ii) we use isotropic variance inside the subspace (as opposed to general Gaussian). Also, using our density model it turns out that to optimize Φ with respect to Σ n , it should be of the form 1 Σ n = β -1 n I d . Therefore we work with β n from now on. Using our density model and g ns = g nκ s and x ns = x nµ s we can write [START_REF] Tipping | Mixtures of probabilistic principal component analysers[END_REF] as:

Φ = ns q ns - d 2 log β n -log q ns - e ns 2σ 2 s - v s 2 dβ -1 n + g ns g ns ρ s + 1 (8) -D log σ s + d 2 log v s ρ s + 1 + log p s + const. with e ns = x ns -α -1 s Λ s R s g ns 2 and v s = ρ s + 1 σ 2 s ρ s α 2 s , (9)
where v s is the inverse variance of p(g | s, x) and e ns is the squared distance between x n and g n mapped into the data space by component s.

Optimization: To optimize Φ we use an EM-style algorithm, a simplified version of the algorithm provided in [START_REF] Roweis | Global coordination of local linear models[END_REF]. The simplifications are: (i) the iterative process to solve for the Λ s , A s is no longer needed; an exact update is possible and (ii) the algorithm no longer involves matrix inversions. The same manner of computation is used: in the E-step, we compute the uni-modal distributions Q(s, g | x n ), parameterized by β n , g n and q ns . Let g n s = E p(g|xn,s) [g] denote the expected value of g given x n and s. We use the following identities:

g n s = κ s + R s Λ s x ns α s ρ s /(ρ s + 1), (10)
D ns = v s 2 dβ -1 n + g n -g n s 2 + d 2 [log β n -log v s ]. (11)
The distributions Q can be found by iterating the fixed-point equations:

β n = s q ns v s , g n = β -1 n s q ns v s g n s , q ns = p ns exp -D ns s p ns exp -D ns ,
where we used p ns = p(s | x n ). In the M-step, we update the parameters of the mixture model. Using notation:

C s = n q ns g ns 2 , E s = n q ns e ns , G s = d n q ns β -1 n , (12)
the update equations are:

κ s = n q ns g n n q ns , µ s = n q ns x n n q ns , α s = C s + G s n q ns (g ns R s Λ s x ns ) , ρ s = D(C s + G s ) d(α 2 s E s + G s ) , σ 2 s = E s + ρ -1 s α -2 s [C s + (ρ s + 1)G s ] (D + d) n q ns
, p s = n q ns ns q ns . Note that the above equations require E s which in turn requires Λ s R s via equations ( 9) and (12). To find Λ s R s we have to minimize:

n q ns e ns = n q ns x ns -α -1 s Λ s R s g ns 2 = - n q ns x ns (Λ s R s )g ns + const.
This problem is known as the 'weighted Procrustes rotation' [START_REF] Cox | Multidimensional Scaling[END_REF]. Let

C = [ √ q 1s x 1s • • • √ q ns x ns ][ √ q 1s g 1s • • • √ q ns g ns ] , with SVD: C = ULΓ ,
where the g ns have been padded with zeros to form D-dimensional vectors, then the optimal Λ s R s is given by the first d columns of UΓ .

Experimental Illustration

To demonstrate the method, we captured 40 × 40 pixel gray valued images of a face with a camera. The face has two degrees of freedom, namely looking updown and left-right. We learned a coordinated mixture model with 1000 images. We used a global PCA projection to 22 dimensions, preserving over 70% of the variance in the data set. We used a latent dimensionality of two and 20 mixture components. We initialized the coordinated mixture model by clamping the latent coordinates g n at coordinates found by Isomap [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF] and clamping the β n at small values for the first 50 iterations. The q ns were initialized uniformly random, and updated from the start. The obtained coordinated mixture model was used to map 1000 'test' images. For each test image x n we approximated p(g

| x n ) with a single Gaussian Q n = arg min Q D KL (Q p(g | x n
)) with a certain mean and standard deviation. In Figure 1 we show these means (location of circle) and standard deviations (radius). To illustrate the discovered parametrization further, two examples of linear traversal of the latent space are given.

Conclusions and Discussion

We showed how a special case of the density model used in [START_REF] Roweis | Global coordination of local linear models[END_REF] leads to a more efficient algorithm to coordinate probabilistic local linear descriptions of a data manifold. The M-step can be computed at once, the iterative procedure to find solutions for a Riccati equation is no longer needed. Furthermore, the update equations do not involve matrix inversion anymore. However, still d singular values and vectors of a D × D matrix have to be found. The application of this method to partially supervised data sets is an interesting possibility and a topic of future research. Another important issue, not addressed here, is that often when we collect data from a system with limited degrees of freedom we actually observe sequences of data. If we assume that the system can vary its state only in continuous manner, these sequences should correspond to paths on the manifold of observable data. This fact might be exploited to find a low dimensional embedding of the manifold. In [START_REF] Verbeek | Procrustes Analysis to Coordinate Mixtures of Probabilistic Principal Component Analyzers[END_REF] we report on promising results of experiments where we used this model to map omni-directional camera images, recorded through an office, to a 2d latent space (the location in the office), where the data was 'supervised' in the sense that 2d workfloor coordinates are known.

Fig. 1 .

 1 Fig. 1. Latent coordinates and linear trajectories in the latent space.

Once we realize that the matrices Vs in[START_REF] Roweis | Global coordination of local linear models[END_REF] are of the form cI d with our density model, it can be seen easily by setting ∂Φ/∂Σn = 0 that Σn = β -1 I d .

Acknowledgment: This research is supported by the Technology Foundation STW (project nr. AIF4997) applied science division of NWO and the technology program of the Dutch Ministry of Economic Affairs.


J J Verbeek

N Vlassis

B Kröse

Fast nonlinear dimensionality reduction with topology preserving networks

We present a fast alternative for the Isomap algorithm. A set of quantizers is fit to the data and a neighborhood structure based on the competitive Hebbian rule is imposed on it. This structure is used to obtain low-dimensional description of the data by means of computing geodesic distances and multi dimensional scaling. The quantization allows for faster processing of the data. The speed-up as compared to Isomap is roughly quadratic in the ratio between the number of quantizers and the number of data points. The quantizers and neighborhood structure are use to map the data to the low dimensional space.

Introduction

In this work we consider high dimensional data X n = {x 1 , . . . , x n } with x i ∈ IR D that lie on a d D dimensional, possibly nonlinear, manifold (plus some noise outside the manifold). The goal is to express the data in the intrinsic coordinates of the manifold. If we are able to map the data to such a low dimensional description then further processing and storage of the data is facilitated, simply because the processing and storage cost is generally at least linear in the dimensionality and also because we may avoid 'the curse of dimensionality'. We combine ideas from Topology Representing Networks [START_REF] Martinetz | Topology representing networks[END_REF], Isomap [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF] and Local Linear Embedding (LLE) [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF].

Our approach first represents the data with a Topology Representing Network (TRN). Then, pairwise distances on the graph between the nodes of the network are computed. The pairwise distances are used to obtain coordinates in IR d of the nodes that mimic the pairwise distances in the graph. These coordinates are obtained by means of Multi Dimensional Scaling (MDS) [START_REF] Cox | Multidimensional Scaling. Number 59 in Monographs on statistics and applied probability[END_REF]. To map the complete data set to IR d local linearity of the manifold is assumed. This allows us to express the data as (convex) linear combinations of nodes of the network in IR D . The same linear combinations are then used to reconstruct the points in IR d .

The virtues of this procedure as compared to Isomap are (i) The computational effort spent by Isomap on constructing the neighborhood graph and the computation of pairwise distances is greatly reduced. The speed up is roughly quadratic in the ratio between the number of quantizers and the number of data points. (ii) In situations where it is unfeasible to store all data, the TRN can be learned online (several methods exist to do so, see for example [START_REF] Fritzke | Handbook of Neural Computation, chapter Unsupervised ontogenetic networks[END_REF]). (iii) Once the TRN is learned, it can be used to map data to and from the latent space. We do not need to learn an explicit mapping, say Radial Basis Function network, from the quantizers to map the other data to the latent space.

Non-linear dimensionality reduction

In this section we briefly discuss three works on nonlinear dimensionality reduction from which ideas are combined in the work presented here.

Topology representing networks Martinetz and Schulten showed [START_REF] Martinetz | Topology representing networks[END_REF] how the simple competitive Hebbian (CHR) rule forms topology representing networks. Let us define Q = {q 1 , . . . , q k } as a set of points, called quantizers, on a manifold M ⊂ IR D . With each quantizer q i a Voronoi set V i is associated in the following manner: V i = {x ∈ IR D : q ix = min j q jx }, where

• denotes the vector norm. The Delaunay triangulation D Q associated with Q is defined as the graph that connects quantizers with adjacent Voronoi sets. 1 The masked Voronoi sets V (M ) i are defined as the intersection of the original Voronoi sets with the manifold M . The Delaunay triangulation D (M ) Q on Q induced by the manifold M is the graph that connects quantizers if the intersection of their masked Voronoi sets is non-empty.

Given a set of quantizers Q and a finite data set X n , the CHR produces a set of edges as follows: (i) For every x i ∈ X n determine the closest and second closest quantizer, respectively q i0 and q i1 . (ii) Include (i 0 , i 1 ) as an edge in E. A set of quantizers Q on M is called dense if for each x on M the triangle formed by x and its closest and second closest quantizer lies completely on M . Obviously, if the distribution of the quantizers over the manifold is homogeneous (the volumes of the associated Voronoi regions are equal), the quantization can be made dense simply by increasing the number of quantizers.

Martinetz and Schulten showed that if Q is dense with respect to M , the CHR produces the induced Delaunay triangulation.

Isomap The Isomap algorithm [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF] finds coordinates in IR d of data that lie on a d dimensional manifold embedded in a D d dimensional space. The aim is to preserve the topological structure of the data, i.e. the Euclidean distances in IR d should correspond to the geodesic distances (distances on the manifold). The algorithm makes use of a neighborhood graph to find the topological structure of the data. The neighborhood graph can be obtained either by connecting all data points that are within some small distance from each other or by connecting each data point to its k nearest neighbors. The algorithm is summarized as follows: (i) Construct neighborhood graph. (ii) Compute the graph distance2 between all data points using a shortest path algorithm, for example Dijkstra's algorithm. (iii) Find low dimensional coordinates by applying MDS on the pairwise distances.

The run time of the Isomap algorithm is dominated by the computation of the neighborhood graph, costing O(n 2 ), and computing the pairwise distances, which costs O(n 2 log n).

Local Linear Embedding The idea underpinning the Local Linear Embedding (LLE) algorithm [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] is the assumption that the manifold is locally linear. It follows that small patches cut out from the manifold in IR D should be approximately equal (up to a rotation, translation and scaling) to small patches on the manifold in IR d . Therefore, local relations among data in IR D that are invariant under rotation, translation and scaling should also be (approximately) valid in IR d . Using this principle, the procedure to find low dimensional coordinates for the data is simple: Express each data point x i as a linear (possibly convex) combination of its k nearest neighbors x i1 , . . . , x i k :

x i = k j=1 w ij x ij + ,
where is the approximation error whose norm is minimized by the weights that are used. Then we want to find coordinates

y i ∈ IR d such that n i=1 y i - k j=1 w ij y ij
2 is minimized. It turns out that the y i can be obtained by finding d eigenvectors of a n × n matrix.

The proposed algorithm

Let us first briefly outline the idea: [START_REF] Cox | Multidimensional Scaling. Number 59 in Monographs on statistics and applied probability[END_REF] The topological structure of the data is captured by a TRN. (2) Low dimensional coordinates for the TRN nodes are computed. [START_REF] Likas | The global k-means clustering algorithm[END_REF] The data is mapped into the low dimensional space using the low dimensional coordinates of the nodes.

(1a) First a set of quantizers that represents the data has to be generated. We propose to use robust variants of the Generalized Lloyd Algorithm (GLA) [START_REF] Likas | The global k-means clustering algorithm[END_REF] to find these. (The GLA is also known as: k-means or LBG algorithm.) (1b) Then, we construct a graph on the quantizer set according to the CHR. Note that the edges emanating from a quantizer can be thought of as indicating the directions in which the manifold extends from this quantizer. (2a) Next, the geodesic distances between the quantizers are estimated by computing graph distances between the quantizers using Dijkstra's algorithm. (2b) We then use MDS to find d dimensional coordinates given these pairwise distances.

Once low dimensional coordinates have been found for the quantizers, we have to assign d dimensional coordinates to the data. We assume that the manifold is locally linear and that the quantizers are distributed densely over the manifold. In this case the manifold is spanned up, locally around a quantizer, by linear combinations of the edges to its neighboring quantizers. Hence, the assumption allows us to write x = j∈N (qx) w j q j + , where q x is the quantizer closest to x, N (q) gives the indices of the neighbors of q (and q itself) in the TRN and is a small error term. Note that the weights w are invariant under scaling, rotation and translation of x and the quantizers in N (q x ). Since the manifold is assumed to be locally linear, the same weights should also describe the data point in the d dimensional space. Therefore, we set the d dimensional coordinates as: y i = j∈N (qi) w i j z j , where z j is the low dimensional coordinate of the j-th quantizer and w i j denotes the weight of quantizer j for reconstructing x i .

To map new data (not present in the training data set) to the low dimensional space, we compute the optimal weights of the closest quantizer and it neighbors in IR D and use the same weights to produce the d dimensional coordinates for the new data point.

Some analysis of the algorithm

Concerning the topology preservation in the nonlinear map provided by MDS we can state a simple result. Given a set of quantizers and a neighborhood graph on them, we define the neighbor-hull of a quantizer q as the convex hull of q and it neighbors in the graph, i.e. {x ∈ IR D : x = j∈N(q) w j q j , w j > 0, j w j = 1}. Two neighbor-hulls are called adjacent if their intersection is non-empty. For the result to hold, we must further constrain the weights w to be non-negative, so that the reconstruction is inside the neighbor-hull, see the comments below.

Lemma 1 Let data points v 1 , v 2 ∈ M ⊂ IR D be touching. If we map these points to their low dimensional representation with our algorithm (using a dense set of quantizers), then they are mapped into adjacent neighbor-hulls in IR d .

Proof Since v 1 , v 2 ∈ M ⊂ IR D are adjacent, they fall in adjacent Voronoi regions, say V 1 and V 2 (associated with quantizers q 1 and q 2 ). Because we use a dense set of quantizers, q 1 , q 2 are connected in the network. Let z 1 , z 2 ∈ IR d be the low dimensional coordinates of q 1 , q 2 . Since the quantizers are adjacent, their neighbor-hulls share the edge (z 1 , z 2 ) in IR d and hence they are adjacent neighbor-hulls.

2 By the same argument we can show that if we map from IR d to IR D using the network obtained by applying the CHR on the data and the quantizers in the low dimensional space, then neighboring points in IR d are mapped into adjacent neighbor-hulls in IR D .

To answer the question of how well the low dimensional representation of the data respected the topological structure of the data, we have two tools. First, as already indicated in [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF], we can measure the correlation between the graph distances for the quantizers and the Euclidean distances in IR d . Second, we can measure the similarity between the networks obtained by applying the CHR in IR d and IR D . In the style of [START_REF] Martinetz | Topology representing networks[END_REF], one might then call the mapping to low dimensionality topology preserving if the two networks are identical.

In order to obtain the weights that reconstruct a point optimally inside the convex hull of the closest quantizer and its neighbors, we have to solve a quadratic programming problem which is computationally demanding. The computational effort can be significantly reduced if also negative weights are allowed, allowing for general linear combinations instead of only convex ones. The computations of these weights is exactly the same as the computation of the weights in LLE. To obtain the optimal weights for x i , up to a re-scaling to satisfy the sum to one constraint, we need to solve Cw = 1. We use w to denote the weight vector. C is a symmetric matrix with C jk = (xq j ) (xq k ), where j, k ∈ N (q i0 ). By adding to C a small multiple of the identity matrix, the reconstructions are drawn toward the center of the convex hull and hence toward the convex hull itself.

To show that the graph distances we use for the MDS mapping converge to manifold distances (how they are interpreted) the same arguments as used for the Isomap algorithm can be used. Note that in particular when the quantizers are assumed to be dense in the sense of Section 2, then is follows directly that the edges have length equal to the manifold distance between the quantizers.

An example on real-world data

We applied our method on data obtained from an omni-directional camera mounted on a Nomad Scout robot. The robot took pictures on 165 positions in an eight by one meter rectangular area in a corridor in our lab. The positions were distributed uniformly over the area, and the orientation of the robot was kept fixed. The omni-directional images were transformed to panoramic images of 400×50 = 20000 pixels. These images we first mapped into a 40 dimensional space by means of Principal Component Analysis, preserving 90% of the total variance in the data. Figure 1 shows the two dimensional coordinates found (by our algorithm and Isomap) for the data together with the quantizer network (connectivity obtained in the 40 dimensional space). The total run time of our algorithm (including the GLA with 30 clusters) was 1.49 seconds, where the Isomap algorithm took 9.62 seconds (five nearest neighbors we used for the neighborhood graph). Both algorithms found an embedding of the data that captures the intrinsic structure in the data and the distance correlation revealed in both cases the correct dimensionality. The speed-up is relatively small in this example, this is due to the small data set used.

Discussion

Advantages (i) If enough data is available and we use c times fewer quantizers as data points we obtain a speed-up, as compared to Isomap3 , of roughly c 2 . (ii) In online learning (or limited memory) settings, we can apply the Dijkstra and MDS steps every once in a while on the online learned TRN, without the Problems (i) In batch-learning settings, the applicability of the method depends on how dense the data-generating manifold is sampled. In very sparsely sampled cases, the network over the quantizers obtained with the CHR will incorrectly represent the manifold. Either we have to many quantizers as compared to data, then the network will consist of many disconnected components. Or we have to few quantizers, which will result in edges connecting areas that are relatively far away on the manifold. (ii) The number of quantizers to use introduces a parameter to be set by the user. However, with the Isomap algorithm one either has to choose the number of neighbors or a maximal distance to build the neighborhood graph, which is unnecessary in the proposed algorithm.

Figure 1 :

 1 Figure 1: The IR 2 coordinates found by our algorithm (left) and Isomap (right).

Two Voronoi sets are called adjacent if their intersection is non-empty.

The graph distance is defined as the minimum length among all paths in the graph that connect the two data points. The length of a path is the sum of the lengths its edges.

In[START_REF] Tenenbaum | Mapping a manifold of perceptual observations[END_REF] also a TRN is used, but a random subset of the data is used for the vertices. Also, the TRN is not used to map the data to IR d .

* This research is supported by the Technology Foundation STW (project nr. AIF4997) applied science division of NWO and the technology program of the Dutch Ministry of Economic Affairs.