With the advent of the big data era, the amount and variety of data are increasing rapidly, making it difficult to estimate the true function. When trying to determine which variables influence a specific value of interest and to what extent, one of the representative statistical methods is the Lasso (least absolute shrinkage and selection operator) regression model (Tibshirani, 1996). This model adds a penalty term to the residual sum of squares in a traditional regression model, facilitating variable selection. In a truncated power spline regression model, if an ℓ₁-norm penalty is imposed on the parameter corresponding to the basis, it can be interpreted as a knot selection problem. When data is given, the spline regression model is important in selecting a significant knot among knot candidates within the range of the predictor, just as variables that influence the response variable are selected through statistical techniques considering all variables. Therefore, this paper focuses on knot selection rather than variable selection, adapting the Lasso method to the context of knot selection and enabling the identification of significant knots.

Much research has been conducted on spline regression models for knot selection (Osborne et al. (1998), Jhong et al. (2017), Jhong and Koo (2019), Oh and Jhong (2021), Lee and Jhong (2021), Lee and Jhong (2022), Lee et al. (2022)). In particular, Osborne et al. (1998), Lee and Jhong (2021), and Lee and Jhong (2022) applied a Lasso penalty and used a truncated power spline as a spline.

Among the various methods to address Lasso regression problems, the alternating direction method of multipliers (ADMM) (Boyd et al., 2011), the coordinate descent algorithm (CDA) (Wright, 2015), and quadratic programming (QP) (Goldfarb and Idnani, 1983) are widely recognized for their efficiency. In this study, we apply the methods to build a Lasso regression model based on truncated power splines and identify which conditions and algorithms are the most efficient in terms of running time and suitability. The methods being compared include ADMM, CDA, and QP. The ADMM algorithm is particularly effective in handling large-scale Lasso regression problems. The CDA has proven to be very efficient for managing high-dimensional data. On the other hand, the QP approach, as detailed in the work of Osborne et al. (2000), leverages the structure of quadratic problems to provide precise solutions for Lasso regression. Despite the availability of these advanced methods, there remains a need for comprehensive comparisons to guide the selection of the most appropriate algorithm for specific situations. By evaluating their performance in terms of suitability and running time, we aim to provide recommendations for their use in different contexts. The results of this comparison will highlight the strengths and weaknesses of each algorithm and offer valuable insights for future research in Lasso regression modeling.

The ADMM algorithm and CDA are directly implemented using the statistical program R, while the iterative algorithm for QP is coded using the solve.QP function in the quadprog package of Weingessel (2011). Additionally, the study utilizes the glmnet function from the R package glmnet, which employs CDA, altering various simulation conditions under the same settings for the four functions. We also examine the fit of truncated power spline regression with a Lasso penalty using two real datasets.

This paper is organized as follows: In Section 2, we present the Lasso regression splines model. The algorithms for comparison are introduced in Section 3, including the glmnet package, which uses the CDA method. In Sections 4 and 5, we validate the performance of our proposed model using simulations and real data. Finally, we summarize the conclusions of this study in Section 6. All numerical analyses were conducted using the statistical software R.

For i = 1, …, n, consider the multiple regression model

where {ɛ_i}′s are independent and are from identical distribution. Here, y_i is i^th sample’s response and x_{i j} is i^th sample’s j^th predictor. Also, δ_j is the parameter related to j^th predictor. Note that for convenience, assuming that y_i is centered, the estimate of the intercept term is 0, so it is omitted. When building a linear regression model using high-dimensional data with many variables, there is a high possibility of numerical instability, making it difficult to find a unique solution for the regression coefficient. In these cases, the Lasso penalty is generally considered. Lasso penalty can be used to build sparse models, which can increase interpretability. The Lasso estimator can be expressed by constraints or can also be expressed by the Lagrange multiplier method. Let δ̂^lasso is (δ^1lasso,…,δ^plasso), then the former is

where t > 0 and the latter is

where λ > 0. As the tuning parameter t decrease and λ increase, the estimated regression coefficient shrinks to 0.

Meanwhile, when given paired data {(xi,yi)}i=1n, we would like to estimate a function f representing the relationship between x_i and y_i. Consider a nonparametric regression model

where x = (x₁, …, x_n) ∈ [0, 1]ⁿ are predictors, y = (y₁, …, y_n) ∈ ℝ ⁿ are responses, and ɛ₁, …, ɛ_n are independent errors with mean zero and a positive variance.

To estimate f, we use truncated power spline basis with K knots t₀ = 0 < t₁ < ··· < t_K < t_K+1 = 1

where α = (α₀, α₁, …, α_d) and β = (β₁, …, β_K). If d = 0,

For d > 0,

Note that (a)₊ is a function that assumes a real value a if a > 0, and 0 otherwise. The number of the initial knots K for the splines is selected and then placed at the equal interval within (0, 1).

Model (2.4) is a form of basis expansion using a truncated power spline basis when a univariate x is given. If these bases are considered as variables, model (2.4) is equivalent to a multiple regression model (2.1). Just as we can select significant variables by applying the Lasso penalty, we select “significant knot(s)” by applying the Lasso penalty in equation (2.5).

The empirical risk of (α, β) is defined as

The objective function to be minimized is given by

for λ > 0 and ||·||₁ is ℓ₁-norm. Let

and then our goal is to estimate lasso regression model f based a truncated power spline.

In equation (2.3), if δ̂_j = 0, the corresponding predictor x_j is not considered in the model. In other words, it can be examined from a variable selection perspective. We can select knots by applying a Lasso penalty to the regression coefficients β_k’s associated with each knot in equation (2.5). Similarly, in equation (2.6), if β̂_k = 0, the corresponding interior knot is not selected. Therefore, just as the simplicity of the model is maintained in the multiple regression model (2.1), the smoothness of the function f, which is the target of estimation, can be adjusted in equation (2.4).

In this study, we apply alternating direction method of multipliers (ADMM), coordinate descent algorithm (CDA), and quadratic programming (QP) to solve the problem (2.6).

3.1. ADMM

The ADMM, which is a relatively newer methodology compared to CDA and QP methodologies, is an algorithm that is intended to blend the decomposability of dual ascent with the superior convergence properties of the method of multipliers. To apply ADMM, the optimization problem (2.6) can be expressed in vector-matrix form using new constraints as follows.

where ||·||₂ is ℓ₂-norm, γ = (γ₁, …, γ_K) ∈ ℝ^K,

and

X is a matrix related to polynomial terms, and P is a matrix generated by a truncated power spline. Then the augmented Lagrangian is

where ξ > 0.

Then, ADMM steps are as follows. For m = 1, …, M,

where M is the number of max iteration. Also, α̃^(m), β̃^(m), γ̃^(m), and ũ^(m) means the m^th updated values of α, β, γ, and u. The vector u ∈ ℝ^K is an additional parameter required when updating.

Herein, we stop iteration when

where ε is a sufficiently small positive number. Additionally, ST(·) is a soft-thresholding operator (Friedman et al., 2010) and is defined as follows:

For λ > 0 and a ∈ ℝ,

Meanwhile, the convergence rate of the algorithm depends on the value of the penalty parameter ξ (Boyd et al., 2011; Xu et al., 2017; McCann and Wohlberg, 2024). If the penalty parameter value is too large, the proportion of minimizing the objective function will be reduced, and if it is too small, the constraints will not be well satisfied.

In our case, we set the following update steps.

where r^(m) is the primal residual at m^th iteration and s^(m) is the dual residual at m^th iteration. Note that α̃⁽⁰⁾, β̃⁽⁰⁾, γ̃⁽⁰⁾ and ũ⁽⁰⁾ are vectors with an element value of 0, which are the initial values of α, β, γ and u, respectively. Also, the initial ξ, ξ⁽⁰⁾ is set λ.

3.2. CDA

The CDA optimizes the objective function with respect to a single coefficient at a time, iteratively cycling through all coefficients until convergence is reached.

For j = 0, 1, …, d and k = 1, 2, …, K, let α̃_j and β̃_k denote the initial/updated values. Specifically,

For notation convenience, we define

Thus, the coordinate descent update is written in the following forms:

Finally, α and β are updated as follows.

It is noteworthy that xi0=1. Here, ri,α(j) and ri,β(k) are the partial residual for α and β, respectively, and

If the difference between the objective functions of the current iteration and the previous iteration is less than _, then the iteration stops.

3.3. QP

The QP is one of the methods for solving mathematical optimization problems involving quadratic functions. Unlike the ADMM algorithm and CDA, it has the characteristic of being able to obtain a solution in single step without repeating the algorithm.

For c ∈ ℝ,

Then, we can represent that c = c₊ − c₋ and |c| = c₊ + c₋. By reparameterizing β_k,

Therefore, the optimization problem is equivalent to minimizing

subject to

in addition to 2 × K non-negativity constraints

Then,

and

where β₊ = {(β_k)₊} and β₋ = {(β_k)₋}. Thus,

Therefore, the minimization problem is equivalent to

subject to

Note that in the QP method, unlike the ADMM algorithm and CDA, when estimating the parameters α and β, the optimization problem (2.6) is transformed into a problem that takes constraints into account, as shown in (2.2). Specifically, the constraints of the optimization problem (3.1) are as follows:

The first constraint corresponds to the Lasso penalty, and the second and third constraints correspond to non-negativity.

To solve the optimization problem (3.1), the solve.QP function of the quadprog package is used in the statistical program R. Dmat, dvec, bvec and Amat are required as arguments for this function. This function solves the optimization problem (− dvec^⊤x + 1/2x^⊤Dmatx) with respect to x having constraints Amat^⊤x ≥ bvec. Therefore, according to the arguments, to solve the optimization problem (3.1), our optimization equation can be expressed as follows:

and

Finally, the solution of b is obtained as

and the estimator is expressed as

3.4. Algorithm details

We implement the algorithms using the R program. To efficiently select tuning parameters, find λ_max, which is the smallest λ among λ’s that assumes all β̂_k to be 0. This is because as λ increases, β̂_k shrinks to 0, and if it increases even a little more than λ_max, all β̂_k = 0 are the same. In our case, the formula for calculating λ_max is as follows (O’Brien, 2016).

where ⟨x,y⟩=Σi=1nxiyi. Thereafter, we generate the λ candidates. After setting λ_ratio, a sufficiently small number, and n_λ, the number of λ candidates, it is split into equal intervals by n_λ from log(λ_max) to log(λ_max × λ_ratio). Next, we transform the candidates back to exponential scale, and replace the last smallest candidate value with 0. When the process of finding λ candidates is completed, we estimate the n_λ models using the ADMM algorithm and CDA.

Here, we specify the maximum number of iterations M and repeat until the convergence criteria are satisfied. After convergence or M iterations, the best model is selected among the n_λ models generated. In statistical data analysis, the criteria for what is “best” exist. Some representative information criteria are, for example, the Akaike information criterion (AIC) (Bozdogan, 1987), the Bayesian information criterion (BIC) (Schwarz, 1978), and cross-validation (CV). The AIC and CV involve the selection of a model with a relatively large variance and small bias compared with the BIC. In general, as the sample size n increases, BIC presents a more sparse model than AIC (Anderson and Burnham, 2004). In addition, BIC can be more effective in selecting the optimal model as the sample size increases because it helps prevent overfitting, and does this better than AIC. Because the purpose of our study is to select a concise model among models that vary depending on tuning parameters, we use BIC, which is defined as follows:

where p_act is the number of non-zero β̂_k’s. We selected the best model that minimized the BIC. Note that Stone (1977) showed that asymptotically minimizing the AIC is equivalent to minimizing the CV value.

Meanwhile, in QP, t is considered as the tuning parameter. Contrary to λ, as t increases, β̂_k approaches least squares estimate (LSE). Therefore, t_max is obtained similarly to λ_max.

Similarly, we determine the number of candidates t, n_t and t_max, and the sufficiently small number t_ratio. Afterwards, it is divided into equal intervals by n_t from log(t_max) to log(t_max × t_ratio), then converted to exponential scale, and returned to the existing scale.

In this section, we analyze the performance using simulated examples based on model (2.1) to compare the four algorithms shown in Section 3: For convenience, we generate the predictor as n sequence and the K knots at regular intervals in the range [0, 1]. ɛ_i is generated from N(0, σ²) for i = 1, …, n. σ² is set differently for each example.

The stopping threshold value ε is utilized to determine the convergence criterion of the algorithm. We set this to 10⁻⁵, to ensure that the algorithm converged sufficiently. M is employed to establish the maximum number of iterations for the algorithm. By setting a limit on the number of iterations, we prevent the algorithm from running indefinitely in cases where it does not converge. The λ_ratio parameter determines the range of λ values by setting the ratio between the maximum and minimum λ values. This adjustment ensures that the regularization parameter is neither excessively large nor too small. The n_λ parameter specifies the number of λ values to be used.

We employ 100 λ values in total to explore a wide range of regularization strengths. We generated models for the various λ values and evaluated their performance based on BIC to select the optimal model. In addition, the simulation is repeated 100 times for each example.

To implement the four algorithms, the statistical program R is used. ADMM and CDA are directly coded, and glmnet used the glmnet function of the glmnet package. In R, glmnet is a package (Friedman et al., 2010; Simon et al., 2013; Friedman et al., 2021; Tay et al., 2023) that fits generalized linear and similar models via penalized maximum likelihood. In this package, glmnet is a function that estimates coefficients using CDA. However, since it is different from the stopping rule considered in ADMM, we also use the glmnet function to compare the performance of the algorithm. Finally, in order to implement the algorithm using QP, parameters are estimated using the solve.QP function of the quadprog package.

Four examples are given as follows:

• Example 1. Constant function

  f(x)=∑hjK(x-xj)K(c)=(1+sgn(c))/2,         c∈ℝsgn(c)={1ifc>00ifc=0-1ifc<0xj=(0.1,0.13,0.15,0.23,0.25,0.40,0.44,0.65,0.76,0.78,0.81)hj=(4,-5,3,-4,5,-4.2,2.1,4.3,-3.1,5.1,-4.2)

• Example 2. Linear function

  f(x)={-3x+3.0if0.0≤x<0.33x+1.2if0.3≤x<0.5-4x+4.7if0.5≤x<0.8x+0.7if0.8≤x≤1.0

• Example 3. Simple cubic function

  f(x)=cos(2πx2)

• Example 4. Complex cubic function

  f(x)=(x(1-x))1/2 sin(2π(1+τ)/(x+τ)),         τ=0.05

Figure 1 represents true functions of Example 1–4, respectively. The top-left plot is for Example 1, which has a constant shape. The top-right plot is for Example 2, which has a continuous fragment linear shape. The bottom-left plot is for Example 3, which has a simple cubic shape using trigonometric function. The bottom-right plot is for Example 4, which has a complex cubic shape. Example 1 and Example 4 utilize functions of Donoho and Johnstone (1994).

To compare the performance of models, we consider the mean squared error (MSE), the mean absolute error (MAE), and the maximum deviation (MXDV) criterion as loss functions. These criteria measure the discrepancy between the true function f and each function f̂. They are expressed as follows.

In addition, we use running time to compare the performance of the models. Comparing running times can help in selecting the most optimal algorithm for each set of different conditions.

In Example 1, n is set to 200 and 400, K to 30 and 40, σ to 0.01 and 0.2, ε to 10⁻⁵, and M to 10, 000. For glmnet, ε is set to 10⁻⁷, and M to 10⁵ by default. Figure 2 shows that the fit of all four algorithms is almost similar to true f. Table 1 shows that the values of MSE, MAE, and MXDV are similar. The fit is similar, however, there is a difference in terms of running time. The running time is long in the order of glmnet < QP < ADMM << CDA. In particular, in the case of CDA, it can be seen that it exceeds 10,000 seconds. Note that the function used in CDA is R user-defined function that implements the CDA update that we derived, whereas glmnet is a function that allows users to easily access Fortran-based code that is highly efficient in numerical calculations in R. That is why glmnet is at least about 86 times faster than CDA in terms of the convergence rate even though it uses the same algorithm. Additionally, CDA user-defined function and glmnet function have a different stopping rule than in ADMM and CDA. For these reasons, it is meaningful to compare whether the two methods are equivalent in terms of fit and the algorithm convergence rate even if the same algorithm is used.

In Example 2, n is set to 200 and 400, K to 20 and 40, σ to 0.05 and 0.15, ε to 10⁻⁵, and M to 10, 000. Looking at Figure 3, the four algorithms appear to fit well, and they seem to estimate the true function better than in Example 1. It is confirmed that the performance indicators of MSE, MAE, and MXDV are all good and similar as seen in Table 2. As in Example 1, the running time is long in the order of glmnet < QP < ADMM << CDA.

In Example 3, n is set to 200 and 400, K to 20 and 40, σ to 0.05 and 0.2, ε to 10⁻⁵, and M to 10, 000. Looking at Figure 4, unlike the previous examples, the fitting performance of QP clearly differs significantly. In Table 3, the MSE value of QP is actually about ten times larger than the MSE values of the other algorithms. The MAE and the MXDV of QP are also about five times larger than the MAE and the MXDV values to the other algorithms. In Figure 4, it can also be observed when x ∈ [0.9, 1.0] that glmnet does not fit well. Therefore, the MSE, MAE, and MXDV values of glmnet are slightly higher compared to those of ADMM and CDA. In terms of running time, similar to the previous examples, CDA takes the longest, followed by ADMM, glmnet, and QP.

In Example 4, n is set to 200 and 400, K to 20 and 40, σ to 0.05 and 0.4, ε to 10⁻⁵, and M to 10, 000. In Figure 5, distinct differences begin to emerge. When x ∈ [0, 0.1], none of the four algorithms performed perfectly. However, when x ∈ [0.1, 0.5], ADMM shows the best fitting result among the four results. In Table 4, the values of MSE, MAE, and MXDV of ADMM are the smallest among the four algorithms. Both CDA and glmnet performed similarly when x ∈ [0.1, 0.5], with their MSE, MAE, and MXDV values also being at similar levels. On the other hand, QP generally does not perform well, with its values of MSE, MAE, and MXDV being the highest among the four algorithms. Moreover, ADMM does not significantly increase the execution time when it changed from a complex cubic function compared to other algorithms.

Meanwhile, in Example 4, ADMM appears to estimate the actual function well. However, CDA, glmnet, and QP show that the fitting results are far from the actual function. This is because, despite using the same basis, ADMM can control the convergence rate with a penalty parameter. In this regard, Figure 6 shows a natural characteristic of using truncated power splines. Figure 6 is a plot expressing the basis created when degrees are 0, 1, 2, and 3, respectively. As the degree increases, the overlap between the basis increases in the plot at large values of support. For this reason, even in the process of estimating the solution, there is correlation between basis, which causes convergence to take a long time or causes problems in estimating the solution.

Note that Figure 7 shows the glmnet fitting results according to the two factors: M and ε for coordinate descent. Figure 7 shows that as M increases and ε decreases, the fitted value gradually converges between x values of 0.4 and 0.6. This observation suggests that achieving convergence in glmnet may require more iterations compared to ADMM under certain conditions.

Conversely, if the penalty parameter is adjusted in ADMM, it is an attractive algorithm that converges faster than CDA. Additionally, for a more objective comparison between ADMM and CDA, we compare the MSE change per iteration for 100 repetitions. In the experiment, the same data as in Example 4 was used, and n = 400, ε = 0.4, and K = 40 were set. Meanwhile, M = 500(m = 1, …, 500) and ε = 0 were set to force the algorithm to run without stopping. Under these settings, in addition to the ADMM and CDA with the ξ update method that we used, we also compare the results when ADMM was fixed to a specific constant value (ξ = 10⁻⁴, 10⁻⁶, 10⁻⁸).

In Figure 8, for each repeat r = 1, …, 100, the x-axis represents the m^th iteration, and the y-axis represents the MSE for f̂ accordingly. The plots were drawn at specific identical λ values for the four ADMM and CDA results. Additionally, the average for 100 lines per method is plotted. First of all, when comparing ADMM(ξ update) and CDA, both algorithms show a fast convergence pattern. Especially, ADMM(ξ update) is absolutely faster than CDA, and when the stopping rule is applied, it actually escapes the iterative algorithm in only 2~3 iterations, whereas CDA shows a monotonically decreasing pattern. Next, to understand the importance of the choice of the penalty parameter, we compare it with the case where ξ is fixed at a specific value instead of updating it as in ADMM. As the ξ value decreases, it shows a pattern of overfitting in the first iteration, and then it gradually converges. Although ADMM(ξ = 10⁻⁴) and ADMM(ξ = 10⁻⁶) are non-monotonous, ADMM(ξ = 10⁻⁴) seems to converge in about the 200th iteration, and ADMM (ξ = 10⁻⁶) is also expected to converge after a few more iterations. However, ADMM(ξ = 10⁻⁸) seems to need much more iterations to converge. In this example, ADMM(ξ update) shows the best performance. In summary, the convergence rate of ADMM varies greatly depending on ξ, which suggests that the ξ value has a great impact on convergence.

In this section, based on the simulation results, we analyze two real data sets. The first is penny, and the second is WWWusage. Both data sets are R built-in data sets.

5.1. Penny data

We utilize the penny dataset included in the locfit package (Loader and Liaw, 2024). The penny dataset contains 90 observations, each containing two variables: year and thickness. The year variable is a discrete variable representing the year when the coin was minted, while the thickness variable is a continuous variable representing the thickness of the coin in millimeters. In Figure 9, the x-axis corresponds to the predictor variable, year, and the y-axis corresponds to the response variable, thickness. For the fitting, K = 30, M = 10⁴ and ε = 10⁻⁵ are set up for all methods, except for glmnet. The glmnet M is set to 10⁵.

The penny’s thickness is not significantly affected by time and shows little variation. The degree of the function for algorithm estimation is zero. As shown in Figure 9, the plot in the top-left has a general trend of increasing thickness from the mid-1940s to the early 1970s. After that, the thickness decreased to around 53(mm), then gradually increased again. The degree of the function for algorithm estimation is zero in fitting penny data. Using a degree of zero is driven by the simplicity of the dataset. It is more efficient for this simple data structure, minimizing overfitting and enhancing model robustness.

The four plots below in Figure 9, show that, except for glmnet, the other algorithms perform at a similarly moderate level. From the mid-1940s to the early 1970s, they fit the changes of the two large data clusters in a simple shape. From the mid-1970s to 1980s, they don’t capture the sharp decline in the data, instead fitting only the average decrease based on the changes in the data clusters. However, glmnet fits well with the general trend of increasing thickness from the mid-1940s to the early 1970s. Additionally, from the mid-1970s to 1980s, it also fits well for the gradual increase following the sharp decline.

5.2. WWWusage data

We utilize the WWWusage dataset, which is included in the base R package. WWWusage is a time series dataset that tracks the number of users connected to the Internet over a specific period. The WWWusage dataset contains a single variable, usage, which represents the number of users connected to the Internet. The dataset consists of 100 observations, recorded at one-minute intervals. In Figure 10, the x-axis corresponds to the predictor variable, minute and the y-axis corresponds to the response variable, number of users.

As shown in Figure 10, the plot in the top-left has a gradual upward trend until the mid-30^th minute. Afterward, the trend decreases until around the 40^th minute. This is followed by a sharp increase until the mid-50^th minute. Subsequently, a decreasing trend is observed until around the 70^th minute. From that point onward, there is a consistent upward trend until the final observation. Overall, compared to the penny dataset, the curve of the WWWusage dataset is much smoother. In the proposed algorithm, min-max normalization is performed on the ‘minute’ variable to avoid errors in numerical calculations, ensuring the predictor variable is within the range [0, 1]. The degree of the function for algorithm estimation is three in fitting WWWusage data. Using a degree of three is driven by the complexity and variability within the dataset. Higher-degree polynomials allow the model to better adapt to the fluctuations and trends. For the fitting, K = 20, M = 10⁴, except for glmnet, and ε = 10⁻⁵ are set up. M of glmnet is 10⁶. The degree of the function for algorithm estimation is three. In the four plots below in Figure 10, ADMM provides a good fit overall, closely following the data across almost all intervals. Both CDA and glmnet only capture the general increasing trend before the sharp decline in the late 50-minute range. However, after 70 minutes, both glmnet and CDA provide good fits for the increasing trend. The QP does not fit the rapid fluctuations from the beginning to the early 50-minute range. Just following the increasing trend only up to around 40 minutes, it fails to capture the sharp increase and shows a decreasing trend until the 70-minute mark. After that, it only followed the general shape of the increasing trend.

We assumed a nonparametric regression equation to estimate the constant, linear, and cubic functions. To solve the problem (2.6), we adopt various methods: ADMM, CDA, and QP. After generating 100 λ’s and fitting the parameters with each algorithm, the best-fit regression model was selected based on the BIC. This process was repeated 100 times, and their performance indicators, such as average MSE, MAE, MXDV, and running time were compared.

For the simpler forms, the constant and linear functions, the fitting performance of the four algorithms was found to be similar. In terms of running time, glmnet perfomed the best. For the simple cubic function, differences in fitting performance were observed, with glmnet and QP not fitting as well as ADMM or CDA. However, in terms of running time, CDA took approximately 3,000 times longer than ADMM. For the complex cubic function, as the number of knots increased, ADMM performed relatively well in terms of both fitting and running time compared to the other algorithms.

Through the simulation, we found that the CDA algorithm has a slow convergence rate in the case of data with complex nonlinear shapes. It was determined that this was due to an issue of overlapping support that occurred when using a truncated power spline as the basis function. In conclusion, we suggest a direction for selecting appropriate algorithms under various conditions. We recommend using glmnet for simpler forms and ADMM for more complex forms.