Regression-clustering for complex real-world data

[0001] The present invention relates to Regression-Clustering (RC) for complex real-world data, and more specifically to a method and system providing an improved RC process based on K-Harmonic Means (KHM) clustering for determining regression functions with reduced residue error from a clustered dataset.

[0002] Complex distributions of computer data inputs are often modeled using a mixture of simpler distributions. Clustering is one of the mathematical tools used to reveal the structure of this mixture. The same is true of data sets with chosen response variables on which a regression analysis can be run. Without separating clusters having very different response properties, the residue error of a regression function is large. Input variable selection could also be misguided to a higher complexity by the mixture.

[0003] In Regression-Clustering, K (>1) regression functions are simultaneously applied to a dataset to guide the clustering into K subsets. Each subset has a simpler distribution for matching to the subsets guiding function. Each function is regressed on its own subset of data thereby resulting in a much smaller residue error. Both the regressions and the clustering optimize a common objective function.

[0004] Two important data mining techniques include regression on data sets with chosen response variables, and clustering on data sets that do not have response information. An RC process is directed at handling the case in between, e.g., data sets that have response variables but do not contain enough information to guarantee high quality learning. The missing part of the response is essential. Missing information is generally caused by insufficiently controlled data collection, due to a lack of means, a lack of understanding or other reasons. For example, sales or marketing data collected on all customers may not have a label on a proper segmentation of the customers. Clustering processes partition a dataset into a finite number of subsets each containing similar data points. Dissimilarity labeled by the index of the partitions provides additional supervision of the K regressions, running in parallel, so that each regression works on a subset of similar data. The K regressions in turn provide the model of dissimilarity for clustering to partition the data. A “linkage” is a common objective function shared between the regressions and the clustering. Neither can be properly done alone without the other.

[0005] Regression-Clustering is not limited to linear regressions, and, when comparing RC between center-based clustering processes, KM (K-Means), KHM (K-Harmonic Means), and EM (Expectation Maximization), the centers are replaced by regression functions. RC refers to a regression-function-centered clustering process. “Clusterwise Linear Regression” uses linear regression and partitioning of the dataset in a process that locally minimizes the total mean square error over all K-regression. Also developed was an incremental version of the process to facilitate adding new observations into the dataset. The Spath process is based on a KM clustering process.

[0006] DeSarbo did research on “Clustered Linear Regression” using the same linear mixing of Gaussian density functions. The number of clusters in the work of Hennig is treated as unknown. Gaffney and Smyth's work is also based on an EM clustering process. Gaffney and Smyth showed applications of regression clustering on video stream data to reveal movements in image sequences.

[0007] Regression-Clustering finds real-world, practical or industrial application in many situations. In economics, demand curves help people to optimize pricing, see Varian, H. R. (1992), “Microeconomic Analysis,” W. W. Norton & Company; 3rd edition. Better understanding of demand curves also helps companies to design multiple models of a product family to fully deploy the area under the demand curves in different segments of a market. Finding the best market segmentation has to be related to the objective that regression is trying to optimize. Regression-Clustering can accomplish both tasks in an integrated process.

[0008] The design of marketing campaigns and offering purchase incentives needs proper segmentation of customers. Without it, marketing campaigns and purchase incentives are blindly given to all potential customers as whole, which is wasteful and less effective. Regression analysis on past marketing campaign data seeks to provide a relationship between an effect and a campaign strategy, e.g., an increase of sales, profit, market share, etc., versus the amount, area or form of the investment, or other. But without proper customer segmentation, regression results are sub-optimal. Regression-Clustering is again a better mathematical tool because Regression-Clustering optimizes both regression and customer-segmentation with a common objective.

[0009] In measuring-device calibrations, regression is run on sampled data to calibrate the device's parameters. However, the accuracy of device may depend on many other factors, some of them may not be controllable or even well understood. The data collected using these devices has missing information, which can be handled by Regression-Clustering. These missing variables can be regarded as either missing input variables or missing response variables. Missing input variables may also be handled by Regression-Clustering in certain situations.

[0010] Many measuring devices work with single-use measuring agents. The manufacturing variations of the measuring agents from different batches are handled by a code, which selects the best set of parameters among multiple sets pre-calibrated and stored in the device. Such code design is based on many runs of regressions on different batches, a costly and time consuming process. Regression-Clustering can optimize both the regression and the clustering (code design) in one step without human intervention, which means significant savings in both time and labor.

[0011] Static or video images can include regions of continuous changes and boundaries of sudden changes in color. A static image can be treated as a mapping from a two-dimensional space to the three-dimensional RGB color-space image:[a,b]×[c,d]→[0,255]×[0,255]×[0,255]. Similarly, a video image can be treated as a mapping from three-dimensional space to another three-dimensional space, video:[a,b]×[c,d]×T→[0,255]×[0,255]×[0,255]. Regression-Clustering is capable of automatically identifying the regions of continuous change and assigning a regression function, which interpolates that part of the image. Both image segmentation and interpolation can be done by Regression-Clustering.

[0012] Previous work on RC used K-Means (KM) and Expectation Maximization (EM) in RC processes, these RC processes have the same well-known problem of being sensitive to the initialization of the regression functions, and the K-Means and EM being sensitive to the initialization of the centers. Previously, a center-based clustering process using K-Harmonic Means has been developed Zhang, B., Hsu, M., Dayal, U. (2000), “K-Harmonic Means”, Intl. Workshop on Temporal, Spatial and Spatio-Temporal Data Mining, Lyon, France Sept. 12; Zhang, B. (2001), “Generalized K-Harmonic Means—Dynamic Weighting of Data in Unsupervised Learning,”, the First SIAM International Conference on Data Mining (SDM'2001), Chicago, USA, Apr. 5-7.

[0013] A KHM center-based clustering process, described in U.S. Pat. No. 6,584,433, issued to the present Assignee, is much less sensitive to initialization of centers than both K-Means and EM. U.S. Pat. No. 6,584,433 describes a harmonic average data clustering method and system. First, a plurality of data points for clustering is received. Next, a number K of clusters is also received. Then, K center points are initialized. For each center point, a new center position is then determined by using a K-Harmonic Means performance function.

[0014] It has been demonstrated through a number of experiments on randomly generated data sets that KHM converges to a better local optimum than K-Means and EM, as measured by a common objective function of K-Means Zhang, B. (2003), “Comparison of the Performance of Center-based Clustering Processes”, the proceedings of PAKDD-03, Seoul, South Korea, April.

[0015] An object of the present invention is to provide an improved method and system for Regression-Clustering (RC) using a K-Harmonic Means (KHM) clustering process. The present invention also provides a new method and system for RC which uses an improved RC KHM process in addition to the possibility of using the existing family of RC processes.

[0016] Briefly, a computer embodiment of the present invention determines regression functions from a data input using K-Harmonic Means (KHM) regression clustering (RC). Processing such data includes (1) selecting K regression functions ƒ₁, . . . , ƒ_K; (2) associating an i-th data point from the dataset with a k-th regression function using a soft membership function; (3) providing a weighting to each data point using a weighting function to determine the data point's participation in calculating a residue error; (4) calculating the residue error between the weighted i-th data point and its associated regression function; and, (5) iterating to minimize any total residue error. The determined regression functions can be applied in data mining, prediction, calibration, segmentation or response analysis.

[0017] An advantage of the present invention is a computer system is provided for RC process based on a K-Harmonic Means clustering process with other existing RC processes based on K-Means and EM.

[0018] Another advantage of the present invention is a computer system is provided for interpretations of the K-regression functions as a predictor and its combination with a K-way classifier.

[0019] These and other objects and advantages of the present invention will no doubt become obvious to those of ordinary skill in the art after having read the following detailed description of the preferred embodiment as illustrated in the drawing figures.

[0031] A particular embodiment of the present invention can be realized using a processing system, an example of which is shown in FIG. 1. In particular, the processing system 100 generally includes at least a processor 102, a memory 104, an input device 106 and an output device 108, connected by a bus 110. An external interface 112 provides for coupling the processing system 100 to a storage device 114 which houses a database 116. The memory 104 can be any form of memory device, for example, volatile or non-volatile memory, solid state storage devices, magnetic devices, etc. The input device 106 can include, for example, a keyboard, pointer device, voice control device, data acquisition card, etc. The output device 108 can include, for example, a display device, monitor, printer, etc.

[0032] In operation, processing system 100 receives a data input 118, processes it, and outputs results in a regression output 120. Such data input is represented in FIGS. 3-12 by the dots in the scatter grams. The regression outputs are represented by one or more curve functions. These can be described by polynomials.

[0033] Such regression functions are computed from the data input using K-Harmonic Means (KHM) regression clustering (RC), as first described by the present inventor in earlier U.S. Patent Applications. Processing such data includes (1) selecting K regression functions ƒ₁, . . . , ƒ_K; (2) associating an i-th data point from the dataset with a k-th regression function using a soft membership function; (3) providing a weighting to each data point using a weighting function to determine the data point's participation in calculating a residue error; (4) calculating the residue error between the weighted i-th data point and its associated regression function; and, (5) iterating to minimize any total residue error.

[0034] FIG. 2 diagrams a K-Harmonic Means regression clustering method embodiment of the present invention for determining regression functions from a computer data input, and is referred to herein by the general reference numeral 200. Such can be executed as software on system 100 (FIG. 1). The process 200 starts with a step 202. A step 204 selects K regression functions ƒ₁, . . . , ƒ_K. At step 206 an i-th data point from the dataset 208 is associated with a k-th regression function obtained at selection step 206 using a soft membership function. At step 210 a weighting is provided to each data point using a weighting function to determine the data point's participation in calculating a residue error. At step 212 the residue error between the weighted i-th data point and its associated regression function is calculated. At step 214 a decision is made whether the total residue error is sufficiently small. If the total residue error is sufficiently small the method 200 can progress to end step 216. If the total residue error is not sufficiently small the steps 204 to 214 are iteratively done to reduce the total residue error until a satisfactory value is achieved.

[0035] Such method is typically embodied as computer software. RC can be built on top of existing regression program libraries and call an existing regression program as a subroutine. The present invention can be applied in complex distributions of real-world data, for example, as a data mining or prediction tool in economics, marketing campaigns, device calibrations, visual image segmentation, etc.

[0036] The processing system 100 can be used to implement an improved Regression-Clustering (RC) process based on a K-Harmonic Means (KHM) clustering process. Given a dataset with supervising responses, Z=(X,Y)={(x_i,y_i)|i=1, . . . , N}, a family of functions Φ={ƒ} and an loss function e( )≧0, regression solves the following minimization problem, $\begin{array}{cc}{f}^{\mathrm{opt}}=\underset{f\in \Phi }{\arg \min }\sum _{i=1}^{N}e\left(f\left(x_i\right),y_i\right)& \left(1\right)\end{array}$
where, Φ is a function class with certain properties to make the optimization problem well defined, such as all polynomials below a certain degree. Usually, $\Phi =\left\{\sum _{l=1}^m{\beta }_{l}h\left(x,a_l\right)❘{\beta }_l\in R,a_l\in R^n\right\},$
linear expansions of simple parametric functions, such as polynomials of degree up to m, Fourier series of bounded frequency, neural networks, RBF, etc. Also usually, e(ƒ(x),y)=∥ƒ(x)−y∥^p, with p=1,2 is most widely used Friedman, J., Hastie, T., and Tibshirani. R. (1998), Additive logistic regression: a statistical view of boosting. Technical report, Department of Statistics, Sequoia Hall, Stanford University, July.
Equation (1) is not effective when the dataset contains a mixture very different response characteristics as shown in FIG. 3, it is much better to find the partitions in the data and learn a separate function on each partition of the dataset as shown in FIG. 4.

[0037] FIG. 3 illustrates a data input 300 that has three likely subgroups 302, 304, and 306. These would ordinarily be lumped together and a single regression 308 would be output.

[0038] FIG. 4 illustrates a data input 400 that has a first subgroup 402 that can be regressed to a function 404, a second subgroup 406 that can be regressed to a function 408, and a third subgroup 410 that can be regressed to a function 412. Thus FIG. 4 shows the extraction of far more information from the data input.

[0039] It can be assumed that there are K partitions in the dataset. Determining the right K number is discussed in a clustering context by Tibshirani, R., Walther, G., and Hastie, T. (2000), “Estimating the Number of Clusters in a Dataset via the Gap Statistic”, see http://www-stat.stanford.edu/˜tibs/research.html; Hamerly and Elkan 2002, see http://www-cse.ucsd.edu/˜ghamerly/academic/papers/icml03.pdf. K can also be determined (or bounded) by other aspects of the original problem.

[0040] In RC processes, K regression functions M={ƒ₁, . . . , ƒ_K}⊂Φ are applied to the dataset, each of which finds its own partition Z_k and regresses on that partition Z_k. Both parts of the process—the K regressions and the partitioning of the dataset—optimize a common objective function. The partition of the dataset can be a “soft” partition given by K density functions defined on the dataset.

[0041] Clusterwise Linear Regression is a simple RC process. The K regressions do not have to be linear. RC-KM solves the following optimization problem, $\begin{array}{cc}\underset{\left\{f_k\right\}⋐\Phi ;\left\{Z_k\right\}}{\min }{\mathrm{Perf}}_{\mathrm{RC}-\mathrm{KM}}=\sum _{k=1}^K\sum _{\left(x_i,y_i\right)\in Z_k}e\left(f_k\left(x_i\right),y_i\right),\mathrm{or}\mathrm{alternatively},\langle \left\{f_k^{\mathrm{opt}}\right\},\left\{Z_k\right\}\rangle =\underset{f\in \Phi ,Z=U_{k=1}^K{z}_k}{\arg \min }\sum _{k=1}^K\sum _{\left(x_i,y_i\right)\in Z_k}e\left(f_k\left(x_i\right),y_i\right)& \left(2\right)\end{array}$
where Z=U_k=1^KZ_k (Z_k I Z_k′=Ø, k≠k′). The optimization is over both the K regression functions and the partition. The optimal partition will satisfy, $\begin{array}{cc}{Z}_k=\left\{\left(x,y\right)\in Z❘e\left(f_k^{\mathrm{opt}}\left(x\right),y\right)\le e\left(f_{k^{\prime }}^{\mathrm{opt}}\left(x\right),y\right)\forall k^{\prime }\ne k\right\},& \left(3\right)\end{array}$
which allows one to replace the function in (2) by, $\begin{array}{cc}{\mathrm{Perf}}_{\mathrm{RC}-\mathrm{KM}}\left(Z,{\left\{f_k\right\}}_{k=1}^K\right)=\sum _{i=1}^N\mathrm{MIN}\left\{e\left(f_k\left(x_i\right),y_i\right)❘k=1,\dots ,K\right\}.& \left(4\right)\end{array}$
The RC-KM Process can be defined as a monotone-convergent process to find a local optimum of (2). Such includes:

[0042] Step 1: Pick K functions ƒ₁⁽⁰⁾, . . . , ƒ_K⁽⁰⁾εΦ randomly, or by any heuristics that are believed to give a good start.

[0043] Step 2: Clustering Phase: In the r-th iteration, r=1, 2, . . . , repartition the dataset as $\begin{array}{cc}{Z}_k^{\left(r\right)}=\left\{\left(x,y\right)\in Z❘e\left(f_k^{\left(r-1\right)}\left(x\right),y\right)\le e\left(f_{k^{\prime }}^{\left(r-1\right)}\left(x\right),y\right)\forall k^{\prime }\ne k\right\}.& \left(5\right)\end{array}$
A tie can be resolved randomly among the winners. Intuitively, each data point is associated with the regression function that gives the smallest approximation error on it. Processically, for r>1, a data point in Z_k^(r−1) is moved to Z_k′^(r) if, $a)e\left(f_{k^{\prime }}^{\left(r-1\right)}\left(x\right),y\right)<e\left(f_k^{\left(r-1\right)}\left(x\right),y\right)\mathrm{and}b)e\left(f_{k^{\prime }}^{\left(r-1\right)}\left(x\right),y\right)\le e\left(f_{k^{″}}^{\left(r-1\right)}\left(x\right),y\right)\mathrm{for}\mathrm{all}{k}^{″}\ne k,k^{\prime }.$
Z_k^(r) inherits all the data points in Z_k^(r−1) that are not moved.

[0044] Step 3: Regression Phase: Run any regression optimization process that gives the following $\begin{array}{cc}{f}_k^{\left(r\right)}=\underset{f\in \Phi }{\arg \min }\sum _{\left(x_i,y_i\right)\in Z_k}e\left(f\left(x_i\right),y_i\right)& \left(6\right)\end{array}$
for k=1, . . . , K. The regression process is selected by the nature of the original problem or other criteria. RC adds no additional constraint on its selection.

[0045] Step 4: Stopping Rule: Run Step 2 and Step 3 repeatedly until there is no more data points changing membership of a subset.

[0046] Steps 2 and 3 do not increase the value of the objective function in equation (2). If any data changes membership in Step 2, the objective function is strictly decreased. Therefore, the process will stop in a finite number of iterations. For reference, see Regularization Girosi, F., Jones, M., and Poggio, T. (1995), “Regularization theory and neural network architectures,” Neural Computation, Vol 7, 219-269; Eldeen, L. (1977), “Processes for the ill-conditioned least square problems,”, BIT, 17, p134-145; and Vapnik, N. V. (1998), “Statistical Learning Theory,” Wiley-Interscience, September to prevent over-fitting, and boosting techniques Schapire, R. E. (1999), “Theoretical views of boosting and applications.” In Tenth International Conference on Processic Learning Theory; and Friedman, J., Hastie, T., and Tibshirani. R. (1998), Additive logistic regression: a statistical view of boosting. Technical report, Department of Statistics, Sequoia Hall, Stanford University, July to improve the quality of the converged results of the regression can also be used on each subset independently.

[0047] Variable selections Montgomery, D. C., Peck, E. A., Vining, G. G. (2001), “Introduction to Linear Regression Analysis, 3rd Edition”, John Wiley & Sons; 3rd edition, April for the K regressions can also be done on each subset independently, with the understanding that an increase in the value of the objective function could result.

[0048] Mean Squared Error (MSE) linear regression with a KM process is one component of some embodiments of the present invention. Assuming {overscore (D)} functions h₁(x), . . . , h_{{overscore (D)}}(x) are chosen as the basis, consider the function class $\Phi =\left\{\sum _{l=1}^D{c}_l{h}_l\left(x\right)❘c_l\in R\right\}.$
To simplify the notations, let
{overscore (x)}=(h₁(x), . . . , h_{{overscore (D)}}(x)) and {overscore (X)}=[{overscore (x)}_i]_{N×{overscore (D)}}.
As an example, for the set of two-variable polynomials up to degree 2, the basis functions are h₁(x)=1, h₂(x)=x₁, h₃(x)=x₂, h₄(x)=x₁², h₅(x)=x₁x₂, h₆(x)=x²². This gives, $X=\left[\begin{array}{cc}{x}_{1,1}& x_{1,2}\\ \cdots & \cdots \\ x_{N,1}& x_{N,2}\end{array}\right]\mathrm{and}$ $\stackrel{\_}{X}=[\begin{array}{cccccc}1& x_{1,1}& x_{1,2}& x_{1,1}^2& x_{1,1}{x}_{1,2}& x_{1,2}^2\\ & \cdots & \cdots & \cdots & \cdots & \cdots \\ 1& x_{N,1}& x_{N,2}& x_{N,1}^2& x_{N,1}{x}_{N,2}& x_{N,2}^2\end{array}].$
With the MSE e(ƒ(x),y)=|ƒ(x)−y|², LinReg-KM minimizes the objective function ${\mathrm{Perf}}_{\mathrm{LinReg}·\mathrm{KM}}\left(Z,{\left\{f_k\right\}}_{k=1}^K\right)=\sum _{i=1}^N\mathrm{MIN}\left\{{{\stackrel{\_}{x}}_i*c_k-y_i}^2❘k=1,\dots ,K\right\}.$
From MSE regression theory, the coefficients of the optimal f is c=({overscore (X)}^T*{overscore (X)})⁻¹{overscore (X)}^T*y and ƒ(x)={overscore (x)}*c.

[0049] With row-partition of Z into K subsets Z₁, . . . , Z_K, matrices {overscore (X)} and Y are row-partitioned accordingly, {overscore (X)}→{overscore (X)}₁, . . . , {overscore (X)}_K and Y→Y₁, . . . , Y_K, the coefficients of the optimal function on the k-th subset is (Step 3 of the RC-KM)
c_k=({overscore (X)}_k^T*{overscore (X)}_k)⁻¹{overscore (X)}_k^T*Y_k.   (7)
The matrix of losses used for the comparisons in Step 2 of RC-KM is
E=[e(ƒ_k(x_i),y_i)]_N×K=abs({overscore (X)}*[c₁, . . . , c_K]−[Y, . . . , Y]).   (8)
There is no need to square the components because squaring is monotone.

[0050] The K-Means clustering process is known to be sensitive to the initialization of its centers. The same is true for RC-KM. Convergence to a poor local optimum has been observed frequently.

[0051] The K-Harmonic Means clustering process showed very strong insensitivity to initialization due to its dynamic weighting of the data points and its non-partitioning membership function Zhang, B. (2001), “Generalized K-Harmonic Means—Dynamic Weighting of Data in Unsupervised Learning,”, the First SIAM International Conference on Data Mining (SDM'2001), Chicago, USA, Apr. 5-7; and Zhang, B. (2003), “Comparison of the Performance of Center-based Clustering Processes”, the proceedings of PAKDD-03, Seoul, South Korea, April.

[0052] An improved method for RC based on a new regression clustering process, RC-KHM_p, can out-perform the RC-KM and RC-EM processes.

[0053] The objective function of RC-KHM_p is defined by replacing the MIN( ) function in equation (4) by a harmonic average HA( ) function. The error function is

[0054]
e(ƒ_k(x_i),y_i)=∥ƒ_k(x_i)−y_i∥^p, p≧2, $\begin{array}{cc}\begin{array}{c}{\mathrm{Perf}}_{\mathrm{RC}-{\mathrm{KHM}}_p}\left(Z,M\right)=\sum _{i=1}^N\underset{1\le k\le K}{\mathrm{HA}}\left\{{f_k\left(x_i\right)-y_i}^p\right\}\\ =\sum _{i=1}^N\frac{K}{\sum _{k=1}^N\frac{1}{{f_k\left(x_i\right)-y_i}^p}}\end{array}& \left(9\right)\end{array}$

The method using the new RC-KHM process includes the steps of:

[0055] Step 1: Pick K functions ƒ_i⁽⁰⁾, . . . , ƒ_K⁽⁰⁾∈Φ randomly or using any heuristic technique believed to offer an improved starting point.

[0056] Step 2: Clustering Phase: In the r-th iteration, let $\begin{array}{cc}{d}_{i,k}=f_k^{\left(r-1\right)}\left(x_i\right)-y_i.& \left(10\right)\end{array}$
The hard partition Z=U_k=1^KZ_k, in RC-KM, is replaced by a “soft” membership function: the i-th data point is associated with the k-th regression function with probability $\begin{array}{cc}p\left(Z_k❘z_i\right)=\frac{d_{i,k}^{p+q}}{\sum _{l=1}^K{d}_{i,l}^{p+q}}.& \left(11\right)\end{array}$
the choice of q will put the regression's error function in L^q-space. See equation (13).

[0057] This is more general than the previous K-Harmonic Means clustering process of the earlier work of the inventor, cited herein, which was limited by only being able to address the situation when q=2.

[0058] For simpler notation, p(Z_k|z_i) and a_p(z_i) in equation (12) are not indexed by q.

[0059] Quantities d_i,k,p(Z_k|z_i), and a_p(z_i) should be indexed by the iteration r, which is also dropped.

[0060] In RC-KHM, not all data points fully participate in all iterations like in RC-KM. Each data point's participation is weighted by $\begin{array}{cc}{a}_p\left(z_i\right)=\frac{\sum _{l=1}^K{d}_{i,l}^{p+q}}{\sum _{l=1}^K{d}_{i,l}^p}.& \left(12\right)\end{array}$
a_p(z_i) is small if and only if z_i is close to one of the functions. The weighting function a_p(z_i) changes in each iteration as the regression functions are updated. If all functions drifted away from a point z_i in the last iteration, a_p(z_i) goes up.

[0061] Step 3: Regression Phase: Run any regression optimization process that gives the following $\begin{array}{cc}{f}_k^{\left(r\right)}=\underset{f\in \Phi }{\arg \min }\sum _{i=1}^N{a}_p\left(z_i\right)p\left(Z_k❘z_i\right){f\left(x_i\right)-y_i}^q& \left(13\right)\end{array}$
for k=1, . . . , K.

[0062] Step 4: Since there is no discrete membership change in RC-KHM, the stopping rule is replaced by measuring the changes to its objective function (9), when the change is smaller than a threshold, the iteration is stopped.

[0063] In Linear Regression with K-Harmonic Means Clustering—LinReg-KHM, q=2 is chosen. Writing equation (13) in matrix form, $\begin{array}{cc}\begin{array}{c}{c}_k^{\left(r\right)}=\underset{c}{\arg \min }{\left(\stackrel{\_}{X}*c-Y\right)}^T*\\ \underset{1\le i\le N}{\mathrm{diag}}\left(a_p\left(z_i\right)p\left(Z_k❘z_i\right)\right)*\left(\stackrel{\_}{X}*c-Y\right)\end{array}& \left(14\right)\end{array}$
and its solution is $\begin{array}{cc}\begin{array}{c}{c}_k^{\left(r\right)}={\left({\stackrel{\_}{X}}^T*{\left[{\stackrel{\_}{x}}_i/{d_{i,k}^{p+2}\left(\sum _{l=1}^K\frac{1}{d_{i,l}^p}\right)}^2\right]}_{\mathrm{Nx}\stackrel{\_}{D}}\right)}^{-1}*\\ {\stackrel{\_}{X}}^T*{\left[y_i/{d_{i,k}^{p+2}\left(\sum _{l=1}^K\frac{1}{d_{i,l}^p}\right)}^2\right]}_{\mathrm{Nx}\stackrel{\_}{D}}\end{array}& \left(15\right)\end{array}$
where d_i,k=∥{overscore (x)}_i*c_k^(r−1)−y_i∥. ([α]_{N×{overscore (D)}} is a matrix of size N×{overscore (D)} with entries α being one of three possibilities: row vectors, column vectors or scalars.) The inversion in equation (15) is on a {overscore (D)}×{overscore (D)} matrix.

[0064] The best of the linear mixing of Gaussian EM clustering process is the natural probability interpretation of its linear mixing (superposition). A brief presentation of RC-EM is included for comparing the performance of all three processes. The objective function for RC-EM is defined as $\begin{array}{cc}\begin{array}{c}{\mathrm{Perf}}_{\mathrm{RC}-\mathrm{EM}}\left(Z,M\right)=-\log \{\prod _{i=1}^N\sum _{k=1}^K\frac{p_k}{\sqrt{{\left(2\pi \right)}^d{\sum }_k}}\\ \mathrm{EXP}\left(-\frac{1}{2}\left(f_k\left(x_i\right)-y_i\right)\sum _k^{-1}{\left(f_k\left(x_i\right)-y_i\right)}^T\right\}\end{array}& \left(16\right)\end{array}$
where d=dim(Y). In case d=1, (ƒ_k(x_i)−y_i) is just a real number and Σ_k⁻¹=1/σ_k².

[0065] The RC-EM recursion is given by $\begin{array}{cc}E\text{-}\mathrm{Step}\text{:}p\left(Z_k^{\left(r\right)}❘z_i\right)=\frac{\frac{p_k^{\left(r-1\right)}}{\sqrt{{\sum }_k}}\mathrm{EXP}\left(-\frac{1}{2}\left(f_k^{\left(r-1\right)}\left(x_i\right)-y_i\right)\sum _{r-1,k}^{-1}{\left(f_k^{\left(r-1\right)}\left(x_i\right)-y_i\right)}^T\right)}{\sum _{k=1}^K\frac{p_k^{\left(r-1\right)}}{\sqrt{{\sum }_k}}\mathrm{EXP}\left(-\frac{1}{2}\left(f_k^{\left(r-1\right)}\left(x_i\right)-y_i\right)\sum _{r-1,k}^{-1}{\left(f_k^{\left(r-1\right)}\left(x_i\right)-y_i\right)}^T\right)}& \left(17\right)\\ M\text{-}\mathrm{Step}\text{:}{p}_k^{\left(r\right)}=\frac{1}{N}\sum _{i=1}^{N}p\left(Z_k^{\left(r\right)}❘z_i\right)& \left(18\right)\\ f_k^{\left(r\right)}=\underset{f\in \Phi }{\arg \min }\sum _{i=1}^{N}p\left(Z_k^{\left(r\right)},z_i\right){f\left(x_i\right)-y_i}^2& \left(19\right)\\ \sum _{r,k}=\frac{\sum _{i=1}^{N}p\left(Z_k^{\left(r\right)}❘z_i\right){\left(f_k^{\left(r\right)}\left(x_i\right)-y_i\right)}^T\left(f_k^{\left(r\right)}\left(x_i\right)-y_i\right)}{N*p_k^{\left(r\right)}}& \left(20\right)\end{array}$

[0066] When MSE linear regression is used, equation (19) can be solved and takes the following special form, while all other equations (16)-(18) and (20) remain the same.
c_k^(r)=({overscore (X)}^T*[p(Z_k^(r),z_i){overscore (x)}_i]_{N×{overscore (D)}})⁻¹*{overscore (X)}^T*[p(Z_k^(r),z_i)y_i]_N×1   (21)
Similarity between equation (21) and LinReg-KHM equation (15), or between equation (21) and LinReg-KM equation (7) is observed.

[0067] The computational cost of one iteration of RC has been compared with the cost of single linear regression on the whole dataset without clustering for all three examples LinReg-KM, LinReg-KHM and LinReg-EM presented herein. Such comparison shows the additional computational cost of switching from single function regression to RC. The comparison is only done for a basic version of regressions without input variable selection or boosting.

[0068] The cost of forming {overscore (X)} is common to both RC and single linear regression.

[0069] In single linear regression, the cost of calculating c=({overscore (X)}^T*{overscore (X)})⁻¹{overscore (X)}^T*Y is the sum of,

[0070] A) {overscore (D)}²*N units for forming {overscore (X)}^T*{overscore (X)},

[0071] B) {overscore (D)}²+{overscore (D)}*N units for forming {overscore (X)}^T*Y

[0072] C) β{overscore (D)}³ for solving ({overscore (X)}^T*{overscore (X)})*c={overscore (X)}^T*Y, β is a small constant.

[0073] where {overscore (D)}=m+1 if D=1, or $\stackrel{\_}{D}=\frac{D^{m+1}-1}{D-1}$
for D>1. A “unit” of calculation here results from multiplying two numbers and adding the result to another number.

[0074] First N≧{overscore (D)}, otherwise the regression has infinite solutions. It is assumed that N>>{overscore (D)}, otherwise the potential of over fitting (and over shoot) is high. In any case the dominate term is O({overscore (D)}²*N).

[0075] Let N_k be the size of the kth cluster, the total cost of K regressions is

[0076] A1) $\sum _{k=1}^K{\stackrel{\_}{D}}^2*N_k={\stackrel{\_}{D}}^2*N$
units for all {overscore (X)}_k^T*{overscore (X)}_k, k=1, . . . , K

[0077] B1) K{overscore (D)}²+{overscore (D)}*N units for all {overscore (X)}_k^T*Y_k and

[0078] C1) Kβ{overscore (D)}³ for solving K linear equations, ({overscore (X)}_k^T*{overscore (X)}_k)*c_k={overscore (X)}_k^T*Y_k

[0079] K is very small and it is not ever expected to be large (say >50).

[0080] The repartition cost for LinReg-KM is O({overscore (D)}*N*K) due to the number of error function evaluations and comparisons. Therefore, the cost of each iteration of LinReg-KM is of the same order of complexity as the simple single function regression.

[0081] The Applicant observed a quick convergence at the start in all experiments but some had a long tail.

[0082] The cost of calculating the repartition probabilities in LinReg-KHM and LinReg-EM are of the same order as the repartition cost in LinReg-KM.

[0083] With input variable selection, not all the variables selected for the single function regression need to appear in the selected variables for each subset. Therefore, the dimensionality of the regression problem on each subset may become lower.

[0084] Regression results are most often used for predictions, y=ƒ(x) is taken as a prediction of the response at a new x∉X. With K regression functions returned by RC, K predictions {ƒ_k(x)}_k=1^K on the same input x are obtained, which is interpreted in this section.

[0085] Assuming that dataset X is Independently and Identically Distributed (IID) sampled from a hidden density distribution P( ). Kernel density estimation Silverman, B. W. (1998), “Density Estimation for Statistics and Data Analysis,” Chapman & Hall/CRC on the K X-projections of Z_k={p(Z_k|z)|z=(x,y)∈Z} (for KHM and EM see equations (11) & (17), for KM they are the real subsets) gives $\begin{array}{cc}\hat{P}\left(x❘X_k\right)=\frac{\frac{1}{N}\sum _{i=1}^{N}p\left(Z_k❘z_i\right)H\left(\frac{x_i-x}{h}\right)}{\hat{P}\left(X_k\right)}\mathrm{with}& \left(22\right)\\ \hat{P}\left(X_k\right)=\frac{1}{N}\sum _{i=1}^{N}p\left(Z_k❘z_i\right).& \left(23\right)\end{array}$
This is a more general subset concept than the classical. LinReg-KM produces real classical subsets, LinReg-KHM and LinReg-EM produce generalized subsets.

[0086] H( ) in equation (22) is a symmetric kernel and h the bandwidth. If the density estimation of each subset is added, the kernel density estimation on the whole dataset is obtained as, $\begin{array}{cc}\hat{P}\left(x\right)=\sum _{k=1}^K\hat{P}\left(x❘X_k\right)\hat{P}\left(X_k\right)=\frac{1}{N}\sum _{i=1}^{N}H\left(\frac{x_i-x}{h}\right).& \left(24\right)\end{array}$
Bayes' inversion gives the probability of x belonging to each subset, $\begin{array}{cc}\hat{P}\left(X_k❘x\right)=\frac{\hat{P}\left(x❘X_k\right)\hat{P}\left(X_k\right)}{\hat{P}\left(x\right)}=\frac{\sum _{i=1}^{N}p\left(Z_k❘z_i\right)H\left(\frac{x_i-x}{h}\right)}{\sum _{i=1}^{N}H\left(\frac{x_i-x}{h}\right)}& \left(25\right)\end{array}$
Let ƒ^%(x) be the random variable prediction which equals ƒ_k(x) with probability P(X_k|x), and the expected value of this prediction is estimated by $\begin{array}{cc}E(f^{\%}\left(x\right)❘\left(x\right)\approx \sum _{k=1}^K{f}_k\left(x\right)\hat{P}\left(X_k❘x\right)=\frac{\sum _{i=1}^N\left[\sum _{k=1}^K{f}_k\left(x\right)p\left(Z_k❘z_i\right)\right]H\left(\frac{x_i-x}{h}\right)}{\sum _{i=1}^{N}H\left(\frac{x_i-x}{h}\right)}.& \left(26\right)\end{array}$
A random variable contains more information than its expectation. Therefore, the RC prediction ƒ^%(x)|x, a random variable, gives more information than its expectation E(ƒ^%(x)|x), which is expected to be close to, but usually not equal to, the single function regression on the whole dataset. Instead of giving a single valued prediction with a large uncertainty, ƒ^%(x)|x gives K possible values each with a much smaller uncertainty. The significant part of the uncertainty is described by the probability distribution {P(X_k|x), k=1, . . . , K}.

[0087] Using the concepts and the relationship that the total variance equals the within-cluster variance plus the between-cluster variance R. Duda & P. Hart, Pattern Classification, 2^nd Ed., Wiley-Interscience (2000), the single value prediction has the total variance. The K-value prediction ƒ^%(x)|x breaks that total variance into the within-cluster variance and the between-cluster variance. The between cluster-variance can be reduced or eliminated if any knowledge outside the dataset helps to choose the k when a new input x is given.

[0088] A classifier, k=C(x), can be trained using the labels provided by the clustering phase of the RC process. In case the false classification rate of C is low, which is not true for all data sets, a prediction on x can be ƒ_C(x)(x).

[0089] The Applicant conducted sets of experiments: Set 1 for visualization of RC, and Set 2 for statistical comparisons of LinReg-KM, LinReg-KHM and LinReg-EM.

[0090] The dimensionality of X is 1, so that 2-dimensional visualization can be presented. Linear regression RC is already demonstrated in FIG. 4. Both quadratic and trigonometric regressions are done. Ploy-KHM, see Zhang, B., Hsu, M., Dayal, U. (2000), “K-Harmonic Means”, Intl. Workshop on Temporal, Spatial and Spatio-Temporal Data Mining, Lyon, France Sept. 12, which performs better for one and two dimensional spaces, is used in this section.

[0091] Referring to FIGS. 5 and 6, the parameters N=600, D=1, K=3 are used. FIG. 5 is the result of simple quadratic regression on a whole dataset comprised of subsets 502, 504, and 506. The single result is a function 508. FIG. 6 represents LinReg-KHM form of the present invention. A subset 602 regresses to a function 604, a subset 606 regresses to a function 608, and a subset 610 regresses to a function 612.

[0092] Referring to FIGS. 7 and 8, the parameters N=1200, D=1, K=3 are used. Φ={a₁sin(6πx)+a₂x+a₃|a_i∈R} and the dataset is a mixture of three subsets generated by three functions in Φ with added Gaussian noise. FIG. 7 shows one regression function applied to the whole dataset. FIG. 8 shows three regression functions being used. Each of the regression functions found a very good approximation of the original functions used to generate the dataset.

[0093] In FIG. 7, a local optimum is shown in FIG. 9 for a single regression function and in FIG. 10 for multiple regression functions. Such hints at how a processes can fail to reach a global optimum. Knowing this can be used to manually correct it, e.g., by providing a special initialization after recognizing a suspected result.

[0094] Referring now to FIG. 11, twelve sets of experiments, with D=2, 4, 6, 8 and K=3, 6, 9, were conducted. In each set, 60 data sets with N=50*D*K were generated by randomly picking N points on K randomly generated hyperplanes and then adding Gaussian noise to the y-components. The regression functions are linear, e.g., hyperplanes. For each dataset, a common initialization of the regression functions was used for all three different processes.

[0095] To make direct comparisons of three processes possible, a common performance measure was used, which was chosen to be the LinReg-KM's objective function in equation (2). After LinReg-KHM and LinReg-EM converged, its own performance measure was discarded and the result re-measured by the LinReg-KM's. Doing so is slightly in favor of LinReg-KM. The notations Perf_KHM/KM and Perf_EM/KM were used for these re-measurements.

[0096] Taking advantage of the known partitions of the synthetic data sets, Perf_baseline, was calculated by running regression on each of the K subsets and adding them up, for comparing against the performance of LinReg-KM and LinReg-KHM. Perf_baseline is close to the global optimum.

[0097] FIG. 11 diagrams the accumulative distribution of selected performance ratios;

• • [0098] squares: LinReg-KHM over LinReg-EM;
  • [0099] (*)'s: LinReg-KHM over the baseline;
  • [0100] (+)'s: LinReg-KM over the baseline,
  • [0101] triangles: LinReg-EM over the baseline.
• [0102] m1=mean of the ratios of LinReg-KHM over LinReg-EM,
• [0103] m2=mean of the ratios of LinReg-KHM over the baseline,
• [0104] m3=mean of the ratios of LinReg-KM over the baseline, and
• [0105] m4=mean of the ratios of LinReg-EM over the baseline.

[0106] Each curve has sixty points from the sixty runs of RC, without interpolation. Four curves in each plot, are frequency-estimations of the accumulative distributions in equations (22)-(25), with v-axis horizontal and prob-axis vertical, $\begin{array}{cc}\mathrm{Pr}\left({\mathrm{Perf}}_{\mathrm{KHM}/\mathrm{KM}}/{\mathrm{Perf}}_{\mathrm{EM}/\mathrm{KM}}<v\right),\mathrm{Pr}\left({\mathrm{Perf}}_{\mathrm{KHM}/\mathrm{KM}}/{\mathrm{Perf}}_{\mathrm{baseline}}<v\right)& \left(22\text{-}23\right)\\ \mathrm{Pr}\left({\mathrm{Perf}}_{\mathrm{RC}-\mathrm{KM}}/{\mathrm{Perf}}_{\mathrm{baseline}}<v\right),\mathrm{Pr}\left({\mathrm{Perf}}_{\mathrm{EM}/\mathrm{KM}}/{\mathrm{Perf}}_{\mathrm{baseline}}<v\right)& \left(24\text{-}25\right)\end{array}$

[0107] The plot of equation (22), in squares, shows how often LinReg-KHM done better then LinReg-EM, with equal performance when the ratio is 1.

[0108] The plot of equation (23), in (*)'s, shows how well LinReg-KHM done against the Perf_baseline, which should be very close to the true optimum. When the value is close to 1, a very good approximation of the global optimum was found.

[0109] The plot of equation (24), in (+)'s and equation (25) in triangles shows how well LinReg-KM and LinReg-EM done against the Perf_baseline.

[0110] The x-axis was truncated to make the interesting part of the plot (near 1) more readable.

[0111] In addition to the plotted distributions in equations (22)-(25), the expectation is also given on each plot, $\begin{array}{cc}\mathrm{m1}\approx E\left(\frac{{\mathrm{Perf}}_{\mathrm{KHM}/\mathrm{KM}}}{{\mathrm{Perf}}_{\mathrm{EM}/\mathrm{KM}}}\right),\mathrm{m2}\approx E\left(\frac{{\mathrm{Perf}}_{\mathrm{KHM}/\mathrm{KM}}}{{\mathrm{Perf}}_{\mathrm{baseline}}}\right),\mathrm{m3}\approx E\left(\frac{{\mathrm{Perf}}_{\mathrm{EM}/\mathrm{KM}}}{{\mathrm{Perf}}_{\mathrm{baseline}}}\right).& \left(26\right)\end{array}$

[0112] Except for K=3 and D=2, LinReg-KHM done the best among the three. As K and D increase, the performance gaps become larger.

[0113] LinReg-EM done better than LinReg-KM on average for all K and D. Such is due to the low dimensionality of the Y-space (dim(Y)=1), where the clustering processes are applied.

[0114] In previous comparisons on the performance of center-based clustering processes Zhang, B. (2003), “Comparison of the Performance of Center-based Clustering Processes”, the proceedings of PAKDD-03, Seoul, South Korea, April, K-means done better than EM on average on data sets with dimensionality>1. The higher the dimensionality of the data, the more K-Means out-performs EM.

[0115] Clustering recovers a discrete estimation of the missing part of the responses and provides each regression function with the right subset of data. A new regression clustering process RC-KHM has been presented herein. It is also observed that LinReg-KHM outperforms both LinReg-EM and LinReg-KM.

[0116] In the general form of RC, the regression part of the process is completely general, no requirements are added by the method using the RC process. Such provides an important insight that (a) RC processes work with any type of regression; and (b) RC can be built on top of existing regression libraries and call an existing regression program as a subroutine.

[0117] Two other advantages of using RC are provided. Regression helps with understanding of a dataset by replacing the dataset with an analytical function plus a residue noise. When the noise is small, the function describes the data well. The compact representation of a dataset by a regression function can also be regarded as data compression, with significantly smaller mean residue noise.

[0118] EM's linear mixing of simple distributions has the most natural probability interpretation. To benefit from both the EM's probability model and the KHM process's robust convergence, it is recommended to run RC-KHM first and use its converged results to initialize RC-EM. RC-KHM does not supply the initial values for p_k^(r) and Σ_r,k. To solve this problem, it is recommended to keep the initial function-centers fixed at the RC-KHM's output for a number of iterations to let the probabilities p_k^(r) and Σ_r,k, if a non-trivial covariance matrix is used, to converge under RC-EM before setting the function-centers free.

[0119] Although the present invention has been described in terms of the presently preferred embodiments, it is to be understood that the disclosure is not to be interpreted as limiting. Various alterations and modifications will no doubt become apparent to those skilled in the art after having read the above disclosure. Accordingly, it is intended that the appended claims be interpreted as covering all alterations and modifications as fall within the true spirit and scope of the invention.