Product Design Time Forecasting by Kernel-Based Regression with Gaussian Distribution Weights

Abstract

There exist problems of small samples and heteroscedastic noise in design time forecasts. To solve them, a kernel-based regression with Gaussian distribution weights (GDW-KR) is proposed here. GDW-KR maintains a Gaussian distribution over weight vectors for the regression. It is applied to seek the least informative distribution from those that keep the target value within the confidence interval of the forecast value. GDW-KR inherits the benefits of Gaussian margin machines. By assuming a Gaussian distribution over weight vectors, it could simultaneously offer a point forecast and its confidence interval, thus providing more information about product design time. Our experiments with real examples verify the effectiveness and flexibility of GDW-KR.

1. Introduction

Product design is a complex and dynamic process, and its duration is affected by a number of factors, most of which are of fuzzy, random and uncertain characteristics. As product design tasks occur in different companies, uncertain characteristics may vary from product to product. The heteroscedasticity thus constitutes another important feature of product design. The mapping from the factors to design time is highly nonlinear, and it is impossible to describe this mapping relationship by definite mathematical models. The degree of reasonability of the supposed distribution of product design time is a key factor in product development control and decisions [1,2,3].

The triangular probability distribution was chosen by Cho and Eppinger [1] to represent design task durations, and a process modeling and analysis technique for managing complex design projects was proposed by using advanced simulation. However, if the assumed distribution of design activity durations does not reflect the true state, the proposed algorithm may fail to obtain ideal results. Yan and Wang [2] proposed a time-computing model with its corresponding design activities in concurrent product development process. Yang and Zhang [3] presented an evolution and sensitivity design-structure matrix to reflect overlapping and their impact on the degree of activity sensitivity and evolution in the process model, and the model can be used for better project planning and control by identifying overlapping and risk for process improvements, but with the two algorithms mentioned above, normal duration of each design activity should be determined before the algorithm is executed, and if activity durations are incompatible with the actual ones, the proposed algorithm may fail to function well. Apparently, the accuracy of predetermined design time is crucial to the planning and controlling of product development processes.

Traditionally, approximate design time is analyzed by means of qualitative approaches. With the rapid development of computer and regression techniques, new forecast methods keep emerging. Bashir and Thomson [4] came up with a modified Norden model to estimate project duration in conjunction with the effort-estimation model. Griffin [5] related the length of the product development cycle to project, process and team structure factors by a statistical method, and quantified the impact of project newness and complexity on the increasing length of development cycle, but with no proposal for design time forecasts. Jacome and Lapinskii [6] developed a model to forecast electronic product design efforts based on a structure and process decomposition approach. Only a small portion of the time factors, however, are taken into account by the model. Xu and Yan [7] proposed a design-time forecast model based on a fuzzy neural network, which exhibits good performance when the sample data are sufficient. However, only a small number of design cases are available to a company, which weakens the validity of the fuzzy neural network. Therefore, a novel approach should be adopted.

Recently, kernel methods have been identified as one of the leading means for pattern classification and function approximation, and successfully applied in various fields [8,9,10,11,12,13,14]. Support vector machine (SVM), initially developed by Vapnik for pattern classification, is one of the most used models. With the introduction of the ε-insensitive loss function, SVM has been extended in use to solve nonlinear regression problems, and thus is also called support vector regression (SVR). ε-insensitive loss functions contribute to the sparseness property of SVR, but the value of ε, chosen a priori, is hard to determine. A new parameter v was then introduced and v-SVR proposed, whereby v controls the number of support vectors and training errors [11]. v-SVR has overcome the difficulty of ε determination. In recent years, much research has been done on kernel methods. Kivinen et al. considered online learning in a reproducing kernel Hilbert space in [15]. Liu et al. [16] proved that the kernel least-mean-square algorithm can be well posed in reproducing kernel Hilbert spaces without adding an extra regularization term to penalize solution norms as was suggested by [15]. Chen et al. developed a quantized kernel least mean square algorithm based on a simple online vector quantization method in [17], and proposed the quantized kernel least squares regression in [18]. Wu et al. [19] derived the kernel recursive maximum correntropy in kernel space and under the maximum correntropy. Furthermore, by combining fuzzy theory with v-SVR, Yan and Xu [20] proposed Fv-SVM to forecast the design time, which could be used to solve regression problems with uncertain input variables. However, both Fv-SVM and v-SVR assume that the noise level is uniform throughout the domain, or at least, its functional dependency is known beforehand [21]. It is thus clear that the time forecast of product design based on Fv-SVM is deficient simply due to the heteroscedasticity of product design. For better planning and controlling of product development process, any good forecast method is expected to yield not only highly precise forecast values, but also valid forecast intervals.

In terms of Gaussian margin machines [22], the weight vector of binary classifier maintains a Gaussian distribution, and what should be struck for is the least information distribution that classifieds training samples with a high probability. Gaussian margin machines provide the probability that a sample belongs to a certain class. The idea given by Gaussian margin machines is extend to the regression for the forecast of product design time. Shang and Yan [23] proposed Gaussian margin regression (GMR) on the basis of combining Gaussian margin machines and kernel-based regression. However, GMR assumes that the forecast variances are same, which is inconsistent with the heteroscedasticity that exits in design time forecast. Like Fv-SVM, GMR also fails to provide valid forecast intervals. By combining Gaussian margin machine and extreme learning machine [24,25], a confidence-weighted extreme learning machine was proposed for regression problems of large samples [26].

The present study adopts the kernel-based regression with Gaussian distribution weights (GDW-KR) by combining Gaussian margin machines with the kernel-based regression, aiming to solve problems of small samples and heteroscedastic noise in design time forecasting, providing both forecast values and intervals. Inheriting the merits of Gaussian margin machines, GDW-KR maintains a Gaussian distribution over weight vectors, seeking the least information distribution that will make each target be included in its corresponding confidence interval. The optimization problem of GDW-KR is simplified, and an approximate solution of the simplified problem is obtained by using the results of regularized kernel-based regression. On the basis of this model, a forecast method for product design time and its relevant parameter-determining algorithm are then put forward.

The rest of this paper is organized as follows: Gaussian margin machines are introduced in Section 2. GDW-KR and the method for solving the optimization problem are described in Section 3. In Section 4, the application in injection mold design is presented, and GDW-KR is then compared with other models. An extended application of GDW-KR is also given. Section 5 draws the final conclusions.

2. Gaussian Margin Machines

Suppose the samples ${\left\{({\mathit{x}}_i,y_i)\right\}}_{i=1}^l$, where ${\mathit{x}}_i\in {\mathbf{R}}^m$ is a column vector and $y_i\in \left\{-1,1\right\}$ is a scalar output. The weight vector $\mathit{w}$ of a linear classifier is supposed to follow a multivariable normal distribution $N_m({\mathit{\mu }}_1,{\mathit{\Sigma }}_1)$ with mean ${\mathit{\mu }}_1\in {\mathbf{R}}^m$ and covariance matrix ${\mathit{\Sigma }}_1\in {\mathbf{R}}^{m\times m}$. For the sample ${\mathit{x}}_i$, we get the normal distribution:

$${{\mathit{x}}_i}^{\mathrm{T}}\mathit{w}\sim N({\mathit{x}}_i^{\mathrm{T}}{\mathit{\mu }}_1,{\mathit{x}}_i^{\mathrm{T}}{\mathit{\Sigma }}_1{\mathit{x}}_i).$$

(1)

The linear classifier is designed to properly classify each sample with a high probability, that is:

$$\mathrm{Pr}(y_i{{\mathit{x}}_i}^{\mathrm{T}}\mathit{w}\ge 0)\ge \rho ,$$

(2)

where $\rho \in (0.5,1]$ is the confidence value.

By combining Equations (1) and (2), we get:

$$\mathrm{Pr}(\frac{y_i{{\mathit{x}}_i}^{\mathrm{T}}\mathit{w}-y_i{{\mathit{x}}_i}^{\mathrm{T}}{\mathit{\mu }}_1}{\sqrt{{\mathit{x}}_i^{\mathrm{T}}{\mathit{\Sigma }}_1{\mathit{x}}_i}}\le \frac{-y_i{{\mathit{x}}_i}^{\mathrm{T}}{\mathit{\mu }}_1}{\sqrt{{\mathit{x}}_i^{\mathrm{T}}{\mathit{\Sigma }}_1{\mathit{x}}_i}})\le 1-\rho .$$

(3)

GMM aims to seek the least informative distribution that classifies the training set with high probability, which is achieved by seeking a multivariable normal distribution $N_m({\mathit{\mu }}_1,{\mathit{\Sigma }}_1)$ with minimum Kullback-Leibler divergence with respect to an isotropic distribution $N_m(\mathbf{0},a{\mathit{I}}_m)$. The Kullback-Leibler divergence between $N_m({\mathit{\mu }}_1,{\mathit{\Sigma }}_1)$ and $N_m(\mathbf{0},a{\mathit{I}}_m)$ is denoted by $D_{\mathrm{KL}}\left(N_m({\mathit{\mu }}_1,{\mathit{\Sigma }}_1)\Vert N_m(\mathbf{0},a{\mathit{I}}_m)\right)$ (the subscript KL is the abbreviation of Kullback-Leibler and D is the abbreviation of divergence), and is obtained by calculating:

$$\frac{1}{2}\ln \det (a{\mathit{I}}_m{\mathit{\Sigma }}_1^{-1})+\frac{1}{2}\mathrm{tr}\left({(a{\mathit{I}}_m)}^{-1}({\mathit{\mu }}_1{\mathit{\mu }}_1^{\mathrm{T}}+{\mathit{\Sigma }}_1-a{\mathit{I}}_m)\right).$$

(4)

The optimization problem of GMM is described as:

$$\begin{array}{cc}\underset{{\mathit{\mu }}_1,{\mathit{\Sigma }}_1}{\min }\hfill & D_{\mathrm{KL}}\left(N_m({\mathit{\mu }}_1,{\mathit{\Sigma }}_1)\Vert N_m(\mathbf{0},a{\mathit{I}}_m)\right)\hfill \\ \mathrm{s.t.}\hfill & \mathrm{Pr}(\frac{y_i{{\mathit{x}}_i}^{\mathrm{T}}\mathit{w}-y_i{{\mathit{x}}_i}^{\mathrm{T}}{\mathit{\mu }}_1}{\sqrt{{\mathit{x}}_i^{\mathrm{T}}{\mathit{\Sigma }}_1{\mathit{x}}_i}}\le \frac{-y_i{{\mathit{x}}_i}^{\mathrm{T}}{\mathit{\mu }}_1}{\sqrt{{\mathit{x}}_i^{\mathrm{T}}{\mathit{\Sigma }}_1{\mathit{x}}_i}})\le 1-\rho \hfill \\ & {\mathit{\Sigma }}_1>0,i=1,\cdots ,l.\hfill \end{array}$$

(5)

After omitting the constant terms in the objective function and transforming the constraints of Equation (5), we get:

$$\begin{array}{l}\underset{{\mathit{\mu }}_1,{\mathit{\Sigma }}_1}{\min }\frac{1}{2}\left(-\ln \det {\mathit{\Sigma }}_1+\frac{1}{a}\mathrm{tr}({\mathit{\Sigma }}_1)+\frac{1}{a}{\mathit{\mu }}_1{\mathit{\mu }}_1^{\mathrm{T}}\right)\\ \mathrm{s.t.}{y}_i{\mathit{x}}_i^{\mathrm{T}}{\mathit{\mu }}_1\ge {\Phi }^{-1}(\rho )\sqrt{{\mathit{x}}_i^{\mathrm{T}}{\mathit{\Sigma }}_1{\mathit{x}}_i}\\ {\mathit{\Sigma }}_1>0,i=1,\cdots ,l,\end{array}$$

(6)

where ${\Phi }^{-1}(\rho )$ is the inverse cumulative distribution function of a standard normal distribution. ${\Phi }^{-1}(\rho )$ is further equal to $\sqrt{2}er{f}^{-1}(2\rho -1)$, where $er{f}^{-1}$ denotes the inverse Gauss error function.

Theorem 1. 

The training samples ${\left\{({\mathit{x}}_i,y_i)\right\}}_{i=1}^l$ are given, and a prior distribution over the weight vector $N_m({\mathit{\mu }}_0,{\mathit{\Sigma }}_0)$ is set. Then, for any $\delta \in [0,1]$ and any posterior distribution $N_m({\mathit{\mu }}_1,{\mathit{\Sigma }}_1)$, the following holds with the probability of at least $1-\delta :$

$$\phi \left(N_m({\mathit{\mu }}_1,{\mathit{\Sigma }}_1),D\right)\le C_1\frac{1}{l}{\displaystyle \sum _{i=1}^l\Phi }(-\frac{y_i{\mathit{x}}_i^{\mathrm{T}}{\mathit{\mu }}_1}{\sqrt{{\mathit{x}}_i^{\mathrm{T}}{\mathit{\Sigma }}_1{\mathit{x}}_i}})+C_2\frac{D_{\mathrm{KL}}\left(N_m({\mathit{\mu }}_1,{\mathit{\Sigma }}_1)\Vert N_m({\mathit{\mu }}_0,{\mathit{\Sigma }}_0)\right)+\ln \frac{2l}{\delta }}{l-1},$$

(7)

where $\phi (N_m({\mathit{\mu }}_1,{\mathit{\Sigma }}_1),D)=\text{E}[\phi (\mathit{w},(\mathit{x},y))|(\mathit{x},y)~D,\mathit{w}\sim N_m({\mathit{\mu }}_1,{\mathit{\Sigma }}_1)]$, $\phi (\mathit{w},(\mathit{x},y))$ is $0-1$ loss function, $C_1=1+\sqrt{2}/2$, $C_2=2+\sqrt{2}/2$, and $D$ is the distribution of $(\mathit{x},y)$ [22,27].

Proof of Theorem 1. 

See Appendix. □

3. Kernel-Based Regression with Gaussian Distribution Weights

3.1. Optimization Problem of GDW-KR

A finite number of independent non-duplicate observations ${\left\{({\mathit{x}}_i,t_i)\right\}}_{i=1}^l$ with ${\mathit{x}}_i\in {\mathbf{R}}^m$ and $t_i\in \mathbf{R}$ are considered. A kernel-based regression model approximates the unknown regression function $f(\mathit{x})$ as follows:

$$\widehat{f}(\mathit{x})={\displaystyle \sum _{j=1}^l{w}_{j}k(\mathit{x},{\mathit{x}}_j),}$$

(8)

where $k(\mathit{x},{\mathit{x}}_j)$ is a predefined kernel function, and $\mathit{w}={(w_1,\dots ,w_l)}^{\mathrm{T}}$.

Definition 1. 

(kernel function) A kernel is a function $k$ that for all $\mathit{x}$, $\mathit{z}$ from a space $\chi$ (which needs not be a vector space) satisfies:

$$k(\mathit{x},\mathit{z})=<\varphi (\mathit{x}),\varphi (\mathit{z})>,$$

(9)

where $\varphi$ is a mapping from the space $\chi$ to a Hilbert space $F$ that is usually called the feature space $\varphi :\mathit{x}\in \chi \mapsto \varphi (\mathit{x})\in F$ [28].

By assuming $\mathit{w}\sim N_l(\mathit{\mu },\mathit{\Sigma })$ with $\mathit{\mu }\in {\mathbf{R}}^l$ and the positive definite covariance matrix $\mathit{\Sigma }\in {\mathbf{R}}^{l\times l}$, we maintain a distribution over alternative weight vectors rather than committing to a single specific vector. Let $y_i$ denote the forecasted value by the model for a given observation ${\mathit{x}}_i$, and we obtain:

$$y_i\sim N_l({\mathit{K}}_i\mathit{\mu },{\mathit{K}}_i\mathit{\Sigma }{\mathit{K}}_i^{\mathrm{T}}),$$

(10)

where ${\mathit{K}}_i$ is the ith row of the symmetric kernel matrix $\mathit{K}$, and ${\mathit{K}}_{ij}=k({\mathit{x}}_i,{\mathit{x}}_j)$, $i=1,\dots ,l$, $j=1,\dots ,l$. Weight vectors are required to make the target value be included in the confidence interval of the forecast value. Thus, we have the following constraint conditions:

$$\begin{array}{l}{\mathit{K}}_i\mathit{\mu }-\eta \sqrt{{\mathit{K}}_i\mathit{\Sigma }{\mathit{K}}_i^{\mathrm{T}}}\le t_i,\\ t_i\le {\mathit{K}}_i\mathit{\mu }+\eta \sqrt{{\mathit{K}}_i\mathit{\Sigma }{\mathit{K}}_i^{\mathrm{T}}},i=1,\dots ,l.\end{array}$$

(11)

The confidence interval needs to be large enough to impose a high confidence level. To make the level higher than 95%, $\eta$ should be greater than 1.96 computed by ${\Phi }^{-1}(1-(1-0.95)/2)$. Considering the independence of noise between samples, ${\mathit{K}}_i\mathit{\Sigma }{\mathit{K}}_{j\ne i}^{\mathrm{T}}$ is set to be $0$. Since the row vector ${\mathit{K}}_i$ cannot be a zero vector, we have ${\mathit{K}}_i\mathit{\Sigma }{\mathit{K}}_i^{\mathrm{T}}>0$, where $\mathit{\Sigma }$ is a positive definite matrix. Hence, the covariance matrix of $\mathit{K}\mathit{\Sigma }{\mathit{K}}^{\mathrm{T}}$ should be a positive definite diagonal matrix:

$$\begin{array}{l}{\mathit{K}}_i\mathit{\Sigma }{\mathit{K}}_{j\ne i}^{\mathrm{T}}=0,\\ \mathit{\Sigma }>0,i=1,\cdots ,l,j=1,\cdots ,l,\end{array}$$

(12)

which indicates that kernel matrix $\mathit{K}$ must be invertible because $\text{rank}(\mathit{K}\mathit{\Sigma }{\mathit{K}}^{\mathrm{T}})=l$ and $\text{rank}(\mathit{K}\mathit{\Sigma }{\mathit{K}}^{\mathrm{T}})\le \text{rank}(\mathit{K})\le l$.

Under the constraint conditions (11) and (12), GDW-KR aims at the least informative distribution that has the smallest Kullback-Leibler divergence with respect to an isotropic Gaussian distribution $N_l(\mathbf{0},a{\mathit{I}}_l)$ for some constant scalar $a>0$. Thus, the optimization problem of GDW-KR is expressed as:

$$\begin{array}{l}\underset{\mathit{\mu },\mathit{\Sigma }}{\min }-\frac{1}{2}\ln \det \mathit{\Sigma }+\frac{1}{2a}\mathrm{tr}(\mathit{\Sigma })+\frac{1}{2a}{\mathit{\mu }}^{\mathrm{T}}\mathit{\mu }\\ \mathrm{s.t.}{\mathit{K}}_i\mathit{\mu }-\eta \sqrt{{\mathit{K}}_i\mathit{\Sigma }{\mathit{K}}_i^{\mathrm{T}}}\le t_i,\\ t_i\le {\mathit{K}}_i\mathit{\mu }+\eta \sqrt{{\mathit{K}}_i\mathit{\Sigma }{\mathit{K}}_i^{\mathrm{T}}},\\ {\mathit{K}}_i\mathit{\Sigma }{\mathit{K}}_{j\ne i}^{\mathrm{T}}=0,\\ \mathit{\Sigma }>0,i=1,\cdots ,l,j=1,\cdots ,l.\end{array}$$

(13)

3.2. Simplification of Optimization Problem

In the problem (13), the number of unknown parameters is $l+l(l+1)/2$, which can be lowered by handling properly its constraints. First of all, let us suppose:

$$\mathit{K}\mathit{\Sigma }{\mathit{K}}^{\mathrm{T}}=\mathit{\Lambda },$$

(14)

where $\mathit{\Lambda }=\mathrm{diag}({\lambda }_1^2,\dots ,{\lambda }_l^2)$, and ${\lambda }_i>0$, $i=1,\cdots ,l$. If the diagonal elements of $\mathit{\Lambda }$ are treated as unknown parameters taking the place of $\mathit{\Sigma }$, the number of unknown parameters in the problem (13) is reduced to $2l$. Then, the objective function of Equation (13) is rewritten as:

$$\underset{\mathit{\mu },\mathit{\Lambda }}{\min }-\frac{1}{2}\ln \det ({\mathit{K}}^{-1}\mathit{\Lambda }{\mathit{K}}^{-1})+\frac{1}{2a}{\mathit{\mu }}^{\mathrm{T}}\mathit{\mu }+\frac{1}{2a}\mathrm{tr}({\mathit{K}}^{-1}\mathit{\Lambda }{\mathit{K}}^{-1}).$$

(15)

As $\ln (\det ({\mathit{K}}^{-1}\mathit{\Lambda }{\mathit{K}}^{-1}))=\ln (\det ({\mathit{K}}^{-1}{\mathit{K}}^{-1}\mathit{\Lambda }))$, we have:

$$-\frac{1}{2}\ln \det ({\mathit{K}}^{-1}\mathit{\Lambda }{\mathit{K}}^{-1})=-{\displaystyle \sum _{i=1}^l\ln {\lambda }_i}-\frac{1}{2}\ln \det \mathit{P},$$

(16)

where $\mathit{P}={\mathit{K}}^{-1}{\mathit{K}}^{-1}$. Since $\mathrm{tr}({\mathit{K}}^{-1}\mathit{\Lambda }{\mathit{K}}^{-1})=\mathrm{tr}({\mathit{K}}^{-1}{\mathit{K}}^{-1}\mathit{\Lambda })$ and both ${\mathit{K}}^{-1}$ and $\mathit{K}$ are symmetric and invertible matrices, we obtain:

$$\mathrm{tr}({\mathit{K}}^{-1}{\mathit{K}}^{-1}\mathit{\Lambda })=\frac{1}{2a}{\displaystyle \sum _{i=1}^l{(\mathit{P})}_{ii}}{\lambda }_i^2,$$

(17)

where ${(\mathit{P})}_{ii}>0$. Disregarding the term $-\frac{1}{2}\ln \det \mathit{P}$ in the objective function, problem (13) is rewritten as:

$$\begin{array}{l}\underset{\mathit{\mu },\mathit{\lambda }}{\min }-{\displaystyle \sum _{i=1}^l\ln {\lambda }_i}+\frac{1}{2a}{\displaystyle \sum _{i=1}^l{(\mathit{P})}_{ii}}{\lambda }_i^2+\frac{1}{2a}{\mathit{\mu }}^{\mathrm{T}}\mathit{\mu }\\ \mathrm{s.t.}{\mathit{K}}_i\mathit{\mu }-\eta {\lambda }_i\le t_i,\\ t_i\le {\mathit{K}}_i\mathit{\mu }+\eta {\lambda }_i,\\ {\lambda }_i>0,i=1,\cdots ,l.\end{array}$$

(18)

Assuming ${\lambda }_i=\lambda$ where in $i=1,\cdots ,l$, the problem of GMR is obtained as:

$$\begin{array}{l}\underset{\mathit{\mu },\lambda }{\min }-l\ln \lambda +\frac{1}{2a}{\lambda }^2{\displaystyle \sum _{i=1}^l{\mathit{P}}_{ii}}+\frac{1}{2a}{\mathit{\mu }}^{\mathrm{T}}\mathit{\mu }\\ \mathrm{s.t.}{\mathit{K}}_i\mathit{\mu }-\eta \lambda \le t_i,\\ t_i\le {\mathit{K}}_i\mathit{\mu }+\eta \lambda ,\\ \lambda >0,i=1,\cdots ,l.\end{array}$$

(19)

Comparing the problems (18) and (19) reveals that GMR is a special case of GDW-KR.

3.3. Analysis of Optimization Problem

Proper generalization of GDW-KR can be guaranteed by Theorem 1 based on the two-sided PAC-Bayesian theorem. However, that of GDW-KR is realized here by analyzing Equation (18) based on the empirical Rademacher complexity [29].

Definition 2. 

(empirical Rademacher complexity) Let $G$ be a family of functions mapping from $X$ to $[a,b]$ and $({\mathit{x}}_1,\dots ,{\mathit{x}}_l)$ a fixed sample of size $l$ with elements in $\mathit{x}$. Then, the empirical Rademacher complexity of $G$ with respect to $({\mathit{x}}_1,\dots ,{\mathit{x}}_l)$ is defined as:

$$\widehat{S}(G)=\underset{\sigma }{\mathrm{E}}\left[\underset{g\in G}{\sup }\left|\frac{2}{l}{\displaystyle \sum _{i=1}^l{\sigma }_{i}g({\mathit{x}}_i)}\right|\right],$$

where $\mathit{\sigma }={({\sigma }_1,\dots ,{\sigma }_l)}^{\mathrm{T}}$ with ${\sigma }_i$ s independent uniform random variables taking values in $\{-1,+1\}$ [30].

Theorem 2. 

GDW-KR can be properly generalized, which is guaranteed by keeping the balance between the empirical Rademacher complexity and the fitting error.

Proof. 

The objective function of the problem (18) is rewritten as:

$$-a{\displaystyle \sum _{i=1}^l\ln {\lambda }_i}+\frac{1}{2}{\displaystyle \sum _{i=1}^l{(\mathit{P})}_{ii}}{\lambda }_i^2+\frac{1}{2}{\mathit{\mu }}^{\mathrm{T}}\mathit{\mu }.$$

(20)

Suppose the function set is as follows:

$$Q_c=\{{\displaystyle \sum _{j=1}^l{\mu }_{j}k(\mathit{x},{\mathit{x}}_j)}|\mathit{x}\in {\mathbf{R}}^m,\mathit{\mu }\in {\mathbf{R}}^l,{\mathit{\mu }}^{\mathrm{T}}\mathit{{\rm K}}\mathit{\mu }\le c^2\},$$

(21)

where $c$ is a positive real number. Let $\widehat{S}(Q_c)$ denote the empirical Rademacher complexity of $Q_c$.

Suppose another function set is defined as:

$$H_c=\{<\mathit{\beta },\varphi (\mathit{x})>|\Vert \mathit{\beta }\Vert \le c\},$$

(22)

where $\varphi$ is the feature mapping corresponding to the kernel $k$.

For any $h(\mathit{x})$ in $H_c$, letting $\mathit{\beta }={\displaystyle \sum _{i=1}^l{\mu }_i\varphi ({\mathit{x}}_i)}$ gives:

$$h(\mathit{x})=<\mathit{\beta },\varphi (\mathit{x})>=<{\displaystyle \sum _{i=1}^l{\mu }_i\varphi ({\mathit{x}}_i)},\varphi (\mathit{x})>={\displaystyle \sum _{i=1}^l{\mu }_{i}k(\mathit{x},{\mathit{x}}_i)},$$

(23)

and:

$${\Vert \mathit{\beta }\Vert }^2=<{\displaystyle \sum _{i=1}^l{\mu }_i\varphi ({\mathit{x}}_i)},{\displaystyle \sum _{j=1}^l{\mu }_j\varphi ({\mathit{x}}_j)}>={\displaystyle \sum _{i,j=1}^l{\mu }_i{\mu }_j<\varphi ({\mathit{x}}_i),\varphi ({\mathit{x}}_j)>}={\displaystyle \sum _{i,j=1}^l{\mu }_i{\mu }_{j}k({\mathit{x}}_i,{\mathit{x}}_j)}={\mathit{\mu }}^{\mathrm{T}}\mathit{{\rm K}}\mathit{\mu }.$$

(24)

Then, $H_c$ is a superset of $Q_c$. Based on the derivation in [30], we obtain $\widehat{S}(Q_c)\le \widehat{S}(H_c)$ and the following:

$$\begin{array}{ll}\widehat{S}(H_c)& =\underset{\sigma }{\mathrm{E}}\left[\underset{h\in H_c}{\sup }\left|\frac{2}{l}{\displaystyle \sum _{i=1}^l{\sigma }_{i}h({\mathit{x}}_i)}\right|\right]\hfill \\ & =\underset{\sigma }{\mathrm{E}}\left[\underset{\Vert \mathit{\beta }\Vert \le c}{\sup }\left|\langle \mathit{\beta },\frac{2}{l}{\displaystyle \sum _{i=1}^l{\sigma }_i\varphi ({\mathit{x}}_i)}\rangle \right|\right]\hfill \\ & \le \frac{2c}{l}\underset{\sigma }{\mathrm{E}}\left[\Vert {\displaystyle \sum _{i=1}^l{\sigma }_i\varphi ({\mathit{x}}_i)}\Vert \right]\hfill \\ & =\frac{2c}{l}\underset{\sigma }{\mathrm{E}}\left[{\left(\langle {\displaystyle \sum _{i=1}^l{\sigma }_i\varphi ({\mathit{x}}_i)},{\displaystyle \sum _{j=1}^l{\sigma }_j\varphi ({\mathit{x}}_j)}\rangle \right)}^{1/2}\right]\hfill \\ & \le \frac{2c}{l}{\left(\underset{\sigma }{\mathrm{E}}\left[{\displaystyle \sum _{i,j=1}^l{\sigma }_i{\sigma }_{j}k({\mathit{x}}_i,{\mathit{x}}_j)}\right]\right)}^{1/2}\hfill \\ & =\frac{2c}{l}{\left({\displaystyle \sum _{i=1}^{l}k({\mathit{x}}_i,{\mathit{x}}_i)}\right)}^{1/2}.\hfill \end{array}$$

Then, we have:

$$\widehat{S}(F_c)\le 2c\sqrt{\mathrm{tr}(\mathit{{\rm K}})}/l.$$

(25)

In view of Equation (21), c can be minimized by minimizing ${\mathit{\mu }}^{\mathrm{T}}\mathit{{\rm K}}\mathit{\mu }$. Calculating by Cauchy–Schwarz inequality yields:

$${\mathit{\mu }}^{\mathrm{T}}\mathit{{\rm K}}\mathit{\mu }=\langle \mathit{\mu },\mathit{{\rm K}}\mathit{\mu }\rangle \le \Vert \mathit{\mu }\Vert \cdot \Vert \mathit{{\rm K}}\mathit{\mu }\Vert \le \Vert \mathit{{\rm K}}\Vert \cdot {\Vert \mathit{\mu }\Vert }^2.$$

(26)

Since the kernel function is predefined, $\frac{1}{2}{\mathit{\mu }}^{\mathrm{T}}\mathit{\mu }$ in Equation (20) can reduce the empirical Rademacher complexity of $Q_c$.

Under the constraints of the problem (18), the smaller ${\lambda }_i$, the less the fitting error. The term:

$$-a{\displaystyle \sum _{i=1}^l\ln {\lambda }_i}+\frac{1}{2}{\displaystyle \sum _{i=1}^l{(\mathit{P})}_{ii}}{\lambda }_i^2.$$

prevents ${\lambda }_i$ from getting too small or too large, and thus the model ${\displaystyle \sum _{j=1}^l{\mu }_{j}k({\mathit{x}}^{\ast },{\mathit{x}}_j)}$ is free from overfitting and underfitting the training data. So the term can be taken as a special loss function. Thereby, it can be concluded that proper values of $a$ and $\eta$ guarantee the balance between the empirical Rademacher complexity and the fitting error. Thus, GDW-KR promises a desirable generalization performance. Then, we have Theorem 2. □

Theorem 2 shows that balancing the empirical Rademacher complexity and the fitting loss is consistent with the two-sided PAC-Bayesian theorem for GDW-KR.

3.4. Solution of Optimization Problem

The results of regularized kernel-based regression are used to obtain the approximate solution of the problem (18). Regularized kernel-based regression is described as:

$$\begin{array}{l}\underset{\mathit{\mu },\mathit{\epsilon }}{\min }\frac{1}{2}({\mathit{\mu }}^{\mathrm{T}}\mathit{\mu }+C{\displaystyle \sum _{i=1}^l{\epsilon }_i^2})\\ \mathrm{s.t.}{\mathit{K}}_i\mathit{\mu }-t_i={\epsilon }_i,\\ i=1,\cdots ,l,\end{array}$$

(27)

where $C$ is the regularization parameter.

Let $\overline{\mathit{\mu }}$ be the solution to Equation (27). Using the KKT conditions, $\overline{\mathit{\mu }}$ is analytically computed as:

$$\overline{\mathit{\mu }}=(\frac{\mathit{I}}{C}+{\mathit{{\rm K}}}^{\mathrm{T}}\mathit{{\rm K}}{)}^{-1}{\mathit{{\rm K}}}^{\mathrm{T}}\mathit{t}.$$

(28)

Then, assuming that $\mathit{\mu }$ is known as $\overline{\mathit{\mu }}$ and ignoring the term $\frac{1}{2a}{\mathit{\mu }}^{\mathrm{T}}\mathit{\mu }$ in the objective function, then we rewrite Equation (18) as:

$$\begin{array}{l}\underset{{\lambda }_i}{\min }-\ln {\lambda }_i+\frac{{(\mathit{P})}_{ii}}{2a}{\lambda }_i^2\\ \mathrm{s.t.}{\lambda }_i\ge t_i^{*},\\ {\lambda }_i>0,\end{array}$$

(29)

where $t_i^{*}=\left|{\mathit{K}}_i\overline{\mathit{\mu }}-t_i\right|/\eta$, $i=1,\cdots ,l$. The second derivative of the objective function of the problem (29) is ${{\lambda }_i}^{-2}+{(\mathit{P})}_{ii}/a$ that must be larger than 0 when ${\lambda }_i>0$. Let ${\overline{\lambda }}_i$ be the solution to Equation (29), which is determined by:

$${\overline{\lambda }}_i=\{\begin{array}{ll}\sqrt{a/{(\mathit{P})}_{ii}},\hfill & t_i^{*}\le \sqrt{a/{(\mathit{P})}_{ii}};\hfill \\ t_i^{*},\hfill & t_i^{*}>\sqrt{a/{(\mathit{P})}_{ii}}.\hfill \end{array}$$

(30)

Thus, the algorithm consists of the following steps:

• Step 1: Make independent non-duplicated observations ${\left\{({\mathit{x}}_i,t_i)\right\}}_{i=1}^l$.
• Step 2: Select the kernel function, and choose the proper relevant parameter (s).
• Step 3: Compute ${\mathit{K}}^{-1}$ and $\mathit{P}$.
• Step 4: Solve the problem (27), and let $\overline{\mathit{\mu }}$ be its solution.
• Step 5: Substitute $\mathit{K}$ and $\overline{\mathit{\mu }}$ into Equation (18), and obtain $\overline{\mathit{\lambda }}$ from Equation (30).

For the observation ${\mathit{x}}^{\ast }$, the forecast value is:

$${\mathit{s}}^{\mathrm{T}}\overline{\mathit{\mu }}={\displaystyle \sum _{j=1}^l{\overline{\mu }}_{j}k({\mathit{x}}^{\ast },{\mathit{x}}_j)},$$

(31)

where $\overline{\mathit{\mu }}={({\overline{\mu }}_1,\dots ,{\overline{\mu }}_l)}^{\mathrm{T}}$, and $\mathit{s}$ = ${\left(k({\mathit{x}}^{\ast },{\mathit{x}}_1),\dots ,k({\mathit{x}}^{\ast },{\mathit{x}}_l)\right)}^{\mathrm{T}}$. And, the forecast interval is calculated as:

$$\left[{\mathit{s}}^{\mathrm{T}}\overline{\mathit{\mu }}-{\eta }^{*}\sqrt{{\mathit{s}}^{\mathrm{T}}\overline{\mathit{\Sigma }}\mathit{s}},{\mathit{s}}^{\mathrm{T}}\overline{\mathit{\mu }}+{\eta }^{*}\sqrt{{\mathit{s}}^{\mathrm{T}}\overline{\mathit{\Sigma }}\mathit{s}}\right],$$

(32)

where $\overline{\mathit{\Sigma }}$ = ${\mathit{K}}^{-1}\overline{\mathit{\Lambda }}{\mathit{K}}^{-1}$ and ${\eta }^{*}>0$.

3.5. Kernel Function and Model Selection

The kernel function plays an important role in kernel function methods. There are three common types of kernel functions: linear function, polynomial function and radial basis function (RBF). Many actual applications demonstrate that RBF tends to display its desirable performance under general smoothness assumptions. With no additional knowledge of the data set available, that makes the very reason for our adoption of the kernel function [31]:

$$k(\mathit{x},{\mathit{x}}_j)=\exp \left\{-{\Vert \mathit{x}-{\mathit{x}}_j\Vert }^2/2{\sigma }^2\right\}.$$

(33)

Hyper-parameters also bear heavily on the generalization performance of kernel function methods. Model selection is to seek proper values of hyper-parameters commonly by means of cross-validation and grid search [32]. The $k$-fold cross-validation [12,13] partitions the training data into $k$ disjoint subsets of approximately equal size. A series of $k$ models are then trained, each using a different combination of k − 1 subsets. The model selection criterion, such as the mean squared error, is then evaluated for each model in each case, utilizing the subset of the data not used in training that model. Recently, evolutional algorithms, such as genetic algorithm and particle swarm optimization, have been adopted to guide the parameters selection process [33,34,35,36]. Regularized kernel-based regression uses genetic algorithm to seek the proper values of $\sigma$ and $C$. An individual in genetic algorithm represents a possible parameter combination. The fitness of each individual is calculated by the $k$-fold cross-validation.

4. Experiments

Experiments were performed to verify the effectiveness of the proposed GDW-KR. The models were built using MATLAB 7.7. The quadratic problems involved were solved through the optimization toolbox QP in MATLAB. The experiments were made on a computer with a Win7 32 bit OS running on 3.1-GHz Intel Core i5-3450 with 4 GB RAM.

4.1. Formulation of Product-Design Time Forecast

To validate the proposed method, the design of plastic injection molds is studied. An injection mold is a kind of single-piece-designed product and the design process is usually driven by customer orders. The design process of injection mold is involved in many product development projects. The design time forecast is meaningful for the optimization of the whole product development process.

Factor values of product-design time are obtained by fuzzy measurable house of quality (FM-HOQ) [7]. Suppose that a design order for a kind of injection mold and the specification of the molding product are given to us. Then the customer demands should be analyzed and some useful mold characteristics should be extracted. The technical customer demands are taken into account. Some demands are originally described as quantitative information (e.g., the mold life is 3000 h), while others are expressed as qualitative information (e.g., the molding product precision is high). A unified fuzzy measurement scheme for all these demands is established, five linguistic levels are used [7]. The importance degrees of these demands are also represented by fuzzy weight sets.

For the specific mold design, the designer should specify the grades of membership of demand weights and demand measures, whose assignments can be made based on the customer demands given on the design order, and on the designer’s objective evaluation of the degrees of importance and scope of the demands. A survey-based methodology is applied for identifying engineering characteristics and time factors, which is performed through self-administered questionnaires from several mold companies in Nanjing. Then, nine kinds of engineering characteristics are selected: mold structure, cavity number, wainscot gauge variation, injection pressure, injection capacity, ejector type, runner shape, manufacturing precision and form feature number. Then we can construct a planning FM-HOQ to map and measure characteristics for technical demands. Among the time characteristics with large influencing weights are structure complexity (SC), model difficulty (MD), wainscot gauge variation (WGV), cavity number (CN), mold size (MS) and form feature number (FFN), the first three of which are expressed as linguistic variables and the last three as numerical ones. Here, the influencing weights that indicate the influence degree on product-design time are different from the indexes of importance in FM-HOQs. Figure 1 presents the application procedure of our model.

4.2. Product-Design Time Forecast Based on GDW-KR

In our experiments, 72 sets of molds with corresponding design time were obtained from a typical company. The detailed characteristic data and design time of these molds compose the corresponding patterns, as shown in

Table 1. Training and testing data of injection model design.

Molds		Input Data						Desired Outputs (h)
No.	Name	SC	MD	WGV	CN	MS	FFN
1	Global handle	L	L	L	4	3.1	3	23
2	Water bottle lid	H	L	H	4	0.56	7	45.5
3	Medicine lid	H	M	VL	4	1.5	6	37
4	Footbath basin	VL	VL	VL	1	0.5	3	10
5	Litter basket	L	M	H	1	2.1	12	42.5
6	Plastic silk flower	L	M	M	1	7.1	4	29.5
…	…	…	…	…	…	…	…	…
71	Paper-lead pulley	L	M	H	8	6.1	6	55
72	Winding tray	M	M	VH	1	2.18	7	41.5

. Numerical variables were normalized to be within [0, 1] by:

$${\overline{x}}_i^d=\frac{x_i^d-\min (x_i^d{|}_{i=1}^l)}{\max (x_i^d{|}_{i=1}^l)-\min (x_i^d{|}_{i=1}^l)},$$

(34)

where $l$ denotes the number of samples, $d$ the number of numerical variables, $x_i^d$ the origin value of the $d\mathrm{th}$ number variable, and ${\overline{x}}_i^d$ the normalized value of the $d\mathrm{th}$ number variable. The linguistic variables, VL, L, M, H and VH, were transformed into the crisp values in terms of expertise: 0.1, 0.25, 0.5, 0.75 and 0.95.

First of all, $\eta$ should be determined, mainly based on the confidence level at which the forecast interval includes the target. To make the confidence level higher than 95%, $\eta$ should be greater than 1.96 computed by ${\Phi }^{-1}(1-(1-0.95)/2)$. The value of $\eta$ is then set to 1.96, and the same is true of ${\eta }^{*}$. The target outputs were normalized to be within [0, 1].

The root mean square error (RMSE), the mean absolute percentage error (MAPE) and the mean absolute error (MAE) are three criteria used to optimize model parameters:

$$\begin{array}{c}\mathrm{RMSE}=\sqrt{\frac{1}{l}{\displaystyle \sum _{i=1}^l{(t_i-{\widehat{t}}_i)}^2}},\\ \mathrm{MAPE}=\frac{1}{l}{\displaystyle \sum _{i=1}^l\left|\frac{t_i-{\widehat{t}}_i}{t_i}\right|},\\ \mathrm{MAE}=\frac{1}{l}{\displaystyle \sum _{i=1}^l\left|t_i-{\widehat{t}}_i\right|},\end{array}$$

where ${\widehat{t}}_i$ is the forecast value for ${\mathit{x}}_i$. The underlying assumption for using the RMSE is that the errors are not biased and follow a normal distribution [37]. The MAPE cannot be used if there is a zero value in $\left\{t_1,\mathrm{...},t_l\right\}$, and puts a heavier penalty on negative errors ($t_i<{\widehat{t}}_i$) than on positive errors. The MAE is suitable to be used for uniformly distributed errors. Because model errors are likely to follow a normal distribution rather than a uniform distribution, the RMSE is a better criterion than the MAE [37]. Thus, we apply the RMSE as a criterion for optimizing model parameters.

The whole data set is divided into several subsets. We choose one subset as the testing set and other ones as the training set. The combination of the genetic algorithm and 5-fold cross-validation is implemented to seek its optimal parameters to minimize the RMSE for the training set. In the genetic algorithm, each individual is evaluated by performing 5-fold cross-validation on the training set. After the optimal parameters are obtained, the model is estimated by using the training set. Then, we calculate the forecast values and three criteria for the testing set. This procedure is repeated until each subset has been used once as the testing set. The testing results of the experiments are averaged over disjoint testing sets which cover the entire dataset. The selection ranges of $\sigma$ and $C$ are [0.01, 5] and [0.01, 10⁶] respectively. The value of $a$ was selected from [10⁻⁶, 10⁶].

The whole data set is first divided into six disjoint subsets. When subset 6 is used as the testing set, the optimal combinational parameters of regularized kernel-based regression are selected as $\sigma =$ 2.119 and $C=$ 998746.999, and the optimal parameter of GDW-KR turns out to be $a=$ 910.190. As illustrated by Figure 2, our GDW-KR gives the valid forecast intervals, excluding T1 and T10. In T10, the forecast interval fails to cover its corresponding target value. In T1, the interval range is too large to provide useful information.

Actual forecast values are listed in

Table 2. Forecast results from four different models when using subset 6 as the testing set.

No.	Designed Outputs	Forecast Results
		Fv-SVM	v-SVR	GMR	GDW-KR
T1	31	35.315	31.236	30.134	32.928
T2	41	40.672	39.167	38.186	39.155
T3	62	62.029	63.521	64.075	63.291
T4	34.5	33.313	32.754	30.900	32.232
T5	16	16.424	16.761	16.877	16.156
T6	32.5	32.965	32.418	32.243	32.801
T7	42.5	40.243	38.516	38.566	38.811
T8	16.5	15.394	15.280	15.346	14.708
T9	22	21.521	21.066	19.963	20.417
T10	54.5	46.391	47.324	46.821	47.789
T11	55	54.149	52.771	53.509	54.304
T12	41.5	41.752	39.893	39.883	40.894

for comparison of the models. The RMSE, the MAPE, the MAE and the average testing time are introduced to compare the forecast performance of different models. Here, the testing time means the time that is spent on solving the optimization problem and on obtaining the testing results when the hyper-parameters are given.

Table 3. Error statistics of four forecast models.

Model	Testing Results			Average Testing Time (s)
	RMSE	MAPE	MAE
Fv-SVM	2.374	0.042	1.905	0.781
v-SVR	2.387	0.041	1.814	0.764
GMR	2.549	0.055	2.137	0.572
GDW-KR	2.366	0.041	1.848	0.583

shows the results from four forecast models, which indicate that GDW-KR promises as high precision as other models do, and that GDW-KR can generate the forecast intervals simultaneously, thus facilitating product development to a certain extent.

The whole data set is then divided into 4 disjoint subsets. Figure 3 illustrates the results of GDW-KR from the first 54 training samples, and demonstrates that GDW-KR still performs well.

Table 4. Error statistics of four forecast models from 54 training samples.

Model	Testing Results			Average Testing Time (s)
	RMSE	MAPE	MAE
Fv-SVM	2.156	0.038	1.615	0.770
v-SVR	2.141	0.037	1.617	0.728
GMR	2.205	0.038	1.624	0.568
GDW-KR	2.133	0.037	1.599	0.579

shows error statistics of four forecast models. GDW-KR does provide a satisfactory performance with small samples, and has thus been proved to be of better performance, appropriate to cases with small samples.

4.3. Extended Application of GDW-KR

Besides design time forecast, GDW-KR can also be extended to other regression problems with small samples. The Slump Test dataset, the Machine CPU dataset and the Yacht Hydrodynamics dataset, which are all from the UCI repository [38], are used to evaluate the extended application of GDW-KR. In these datasets, Fv-SVM behaves the same as v-SVR, as there is no fuzzy variable. Thus, the results of Fv-SVM are not presented here. Each dataset is divided into 6 disjoint subsets. In our experiments, both the target output and numerical attributes were normalized to be within [0, 1].

The Concrete Slump Test covers seven input and three output variables as well as 103 data points. The 28-day Compressive Strength is taken as the desired output variable. For the case of the Concrete Slump Test, the results of GDW-KR are compared with those of other two models. Concrete Slump Test results are shown in Figure 4. The three error indices of different three models are given in

Table 5. Error statistics of three forecast models on the Slump Test dataset.

Model	Testing Results			Average Testing Time (s)
	RMSE	MAPE	MAE
v-SVR	0.021	0.044	0.014	0.795
GMR	0.023	0.047	0.015	0.583
GDW-KR	0.019	0.055	0.014	0.607

. On the Slump Test, GDW-KR offers forecast values with high accuracy and forecast intervals with good validity.

For the Machine CPU dataset and the Yacht Hydrodynamics, the error statistics of three forecast models are presented in

Table 6. Error statistics of three forecast models on the Machine CPU.

Model	Testing Results			Average Testing Time (s)
	RMSE	MAPE	MAE
v-SVR	0.038	0.660	0.024	0.941
GMR	0.039	0.638	0.023	0.695
GDW-KR	0.038	0.555	0.023	0.783

and

Table 7. Error statistics of three forecast models on the Yacht Hydrodynamics.

Model	Testing Results			Average Testing Time (s)
	RMSE	MAPE	MAE
v-SVR	0.034	3.585	0.025	1.196
GMR	0.036	3.772	0.027	0.710
GDW-KR	0.034	3.471	0.026	0.894

, respectively. Figure 5 and Figure 6 indicate the forecast results when using subset 6 as the testing set.

5. Conclusions

The control and decision of product development are based on the reasonable degree of the distribution of product design time. In design time forecasting, the problems of small samples and heteroscedastic noise ought to be considered.

This paper has presented a new model of kernel-based regression with Gaussian distribution weights for product-design time forecasts, which combines Gaussian margin machines with kernel-based regression. The kernel method performs well for the problem of small samples. Unlike GMR, which assumes that the covariance matrix of the forecast values in the training set is an identity matrix multiplied by a positive scalar, GDW-KR assumes that this matrix is a positive definite diagonal matrix. GDW-KR is more suitable for addressing the problem of heteroscedastic noise than GMR, and has the advantage of providing both point forecasts and confidence intervals simultaneously.

The plastic injection mold was studied before modeling. For convincing evaluation, experiments with 72 real samples were conducted. Results from them have verified that GDW-KR promises not only as high forecast accuracy as Fv-SVM and v-SVR but forecast intervals crucial to the control and decision of product development. Undoubtedly, GDW-KR benefits from the merits of Gaussian margin machines.