Search

In this paper, we introduce an HPR-LP solver, an implementation of a Halpern Peaceman-Rachford (HPR) method with semi-proximal terms for solving linear programming (LP). The HPR method enjoys the iteration complexity of $O(1/k)$ in terms of the Karush-Kuhn-Tucker residual and the objective error. Based on the complexity results, we design an adaptive strategy of restart and penalty parameter update to improve the efficiency and robustness of the HPR method. We conduct extensive numerical experiments on different LP benchmark datasets using NVIDIA A100-SXM4-80GB GPU in different stopping tolerances. Our solver's Julia version achieves a $\textbf{2.39x}$ to $\textbf{5.70x}$ speedup measured by SGM10 on benchmark datasets with presolve ($\textbf{2.03x}$ to $\textbf{4.06x}$ without presolve) over the award-winning solver PDLP with the tolerance of $10^{-8}$.

In this paper, we study the Aubin property of the Karush-Kuhn-Tucker solution mapping for the nonlinear semidefinite programming (NLSDP) problem at a locally optimal solution. In the literature, it is known that the Aubin property implies the constraint nondegeneracy by Fusek [SIAM J. Optim. 23 (2013), pp. 1041-1061] and the second-order sufficient condition by Ding et al. [SIAM J. Optim. 27 (2017), pp. 67-90]. Based on the Mordukhovich criterion, here we further prove that the strong second-order sufficient condition is also necessary for the Aubin property to hold. Consequently, several equivalent conditions including the strong regularity are established for NLSDP's Aubin property. Together with the recent progress made by Chen et al. on the equivalence between the Aubin property and the strong regularity for nonlinear second-order cone programming [arXiv:2406.13798v1 (2024)], this paper constitutes a significant step forward in characterizing the Aubin property for general non-polyhedral $C^2$-cone reducible constrained optimization problems.

This paper proposes an efficient HOT algorithm for solving the optimal transport (OT) problems with finite supports. We particularly focus on an efficient implementation of the HOT algorithm for the case where the supports are in $\mathbb{R}^2$ with ground distances calculated by $L_2^2$-norm. Specifically, we design a Halpern accelerating algorithm to solve the equivalent reduced model of the discrete OT problem. Moreover, we derive a novel procedure to solve the involved linear systems in the HOT algorithm in linear time complexity. Consequently, we can obtain an $\varepsilon$-approximate solution to the optimal transport problem with $M$ supports in $O(M^{1.5}/\varepsilon)$ flops, which significantly improves the best-known computational complexity. We further propose an efficient procedure to recover an optimal transport plan for the original OT problem based on a solution to the reduced model, thereby overcoming the limitations of the reduced OT model in applications that require the transport map. We implement the HOT algorithm in PyTorch and extensive numerical results show the superior performance of the HOT algorithm compared to existing state-of-the-art algorithms for solving the OT problems.

We introduce ``LAMBDA," a novel open-source, code-free multi-agent data analysis system that that harnesses the power of large models. LAMBDA is designed to address data analysis challenges in complex data-driven applications through the use of innovatively designed data agents that operate iteratively and generatively using natural language. At the core of LAMBDA are two key agent roles: the programmer and the inspector, which are engineered to work together seamlessly. Specifically, the programmer generates code based on the user's instructions and domain-specific knowledge, enhanced by advanced models. Meanwhile, the inspector debugs the code when necessary. To ensure robustness and handle adverse scenarios, LAMBDA features a user interface that allows direct user intervention in the operational loop. Additionally, LAMBDA can flexibly integrate external models and algorithms through our knowledge integration mechanism, catering to the needs of customized data analysis. LAMBDA has demonstrated strong performance on various machine learning datasets. It has the potential to enhance data science practice and analysis paradigm by seamlessly integrating human and artificial intelligence, making it more accessible, effective, and efficient for individuals from diverse backgrounds. The strong performance of LAMBDA in solving data science problems is demonstrated in several case studies, which are presented at \url{https://www.polyu.edu.hk/ama/cmfai/lambda.html}., Comment: 30 pages, 21 figures and 5 tables

Semismooth* Newton methods have been proposed in recent years targeting multi-valued inclusion problems and have been successfully implemented to deal with several concrete generalized equations. In this paper, we show that these executable implementations are exactly the applications of G-semismooth Newton methods for solving nonsmooth equations localized from these generalized equations. This new understanding expands the breadth of G-semismooth Newton methods in theory, and more importantly, facilitates the design and implementation of practical Newton-type algorithms for solving generalized equations.

This paper solves a fundamental open problem in variational analysis on the equivalence between the Aubin property and the strong regularity for nonlinear second-order cone programming (SOCP) at a locally optimal solution. We achieve this by introducing a reduction approach to the Aubin property characterized by the Mordukhovich criterion and a lemma of alternative choices on cones to replace the S-lemma used in Outrata and Ram\'irez [SIAM J. Optim. 21 (2011) 789-823] and Opazo, Outrata, and Ram\'irez [SIAM J. Optim. 27 (2017) 2141-2151], where the same SOCP was considered under the strict complementarity condition except for possibly only one block of constraints. As a byproduct, we also offer a new approach to the well-known result of Dontchev and Rockafellar [SIAM J. Optim. 6 (1996) 1087-1105] on the equivalence of the two concepts in conventional nonlinear programming.

Real-time prediction plays a vital role in various control systems, such as traffic congestion control and wireless channel resource allocation. In these scenarios, the predictor usually needs to track the evolution of the latent statistical patterns in the modern high-dimensional streaming time series continuously and quickly, which presents new challenges for traditional prediction methods. This paper is the first to propose a novel online algorithm (TOPA) based on tensor factorization to predict streaming tensor time series. The proposed algorithm TOPA updates the predictor in a low-complexity online manner to adapt to the time-evolving data. Additionally, an automatically adaptive version of the algorithm (TOPA-AAW) is presented to mitigate the negative impact of stale data. Simulation results demonstrate that our proposed methods achieve prediction accuracy similar to that of conventional offline tensor prediction methods, while being much faster than them during long-term online prediction. Therefore, TOPA-AAW is an effective and efficient solution method for the online prediction of streaming tensor time series.

In this paper, we aim to accelerate a preconditioned alternating direction method of multipliers (pADMM), whose proximal terms are convex quadratic functions, for solving linearly constrained convex optimization problems. To achieve this, we first reformulate the pADMM into a form of proximal point method (PPM) with a positive semidefinite preconditioner which can be degenerate due to the lack of strong convexity of the proximal terms in the pADMM. Then we accelerate the pADMM by accelerating the reformulated degenerate PPM (dPPM). Specifically, we first propose an accelerated dPPM by integrating the Halpern iteration and the fast Krasnosel'ski\u{i}-Mann iteration into it, achieving asymptotic $o(1/k)$ and non-asymptotic $O(1/k)$ convergence rates. Subsequently, building upon the accelerated dPPM, we develop an accelerated pADMM algorithm that exhibits both asymptotic $o(1/k)$ and non-asymptotic $O(1/k)$ nonergodic convergence rates concerning the Karush-Kuhn-Tucker residual and the primal objective function value gap. Preliminary numerical experiments validate the theoretical findings, demonstrating that the accelerated pADMM outperforms the pADMM in solving convex quadratic programming problems.

Unsupervised feature selection has drawn wide attention in the era of big data since it is a primary technique for dimensionality reduction. However, many existing unsupervised feature selection models and solution methods were presented for the purpose of application, and lack of theoretical support, e.g., without convergence analysis. In this paper, we first establish a novel unsupervised feature selection model based on regularized minimization with nonnegative orthogonal constraints, which has advantages of embedding feature selection into the nonnegative spectral clustering and preventing overfitting. An effective inexact augmented Lagrangian multiplier method is proposed to solve our model, which adopts the proximal alternating minimization method to solve subproblem at each iteration. We show that the sequence generated by our method globally converges to a Karush-Kuhn-Tucker point of our model. Extensive numerical experiments on popular datasets demonstrate the stability and robustness of our method. Moreover, comparison results of algorithm performance show that our method outperforms some existing state-of-the-art methods.

Recently, a quaternion tensor product named Qt-product was proposed, and then the singular value decomposition and the rank of a third-order quaternion tensor were given. From a more applicable perspective, we extend the Qt-product and propose a novel multiplication principle for third-order quaternion tensor named gQt-product. With the gQt-product, we introduce a brand-new singular value decomposition for third-order quaternion tensors named gQt-SVD and then define gQt-rank and multi-gQt-rank. We prove that the optimal low-rank approximation of a third-order quaternion tensor exists and some numerical experiments demonstrate the low-rankness of color videos. So, we apply the low-rank quaternion tensor completion to color video inpainting problems and present alternating least-square algorithms to solve the proposed low gQt-rank and multi-gQt-rank quaternion tensor completion models. The convergence analyses of the proposed algorithms are established and some numerical experiments on various color video datasets show the high recovery accuracy and computational efficiency of our methods.

In an era of digital ubiquity, efficient resource management and decision-making are paramount across numerous industries. To this end, we present a comprehensive study on the integration of machine learning (ML) techniques into Huawei Cloud's OptVerse AI Solver, which aims to mitigate the scarcity of real-world mathematical programming instances, and to surpass the capabilities of traditional optimization techniques. We showcase our methods for generating complex SAT and MILP instances utilizing generative models that mirror multifaceted structures of real-world problem. Furthermore, we introduce a training framework leveraging augmentation policies to maintain solvers' utility in dynamic environments. Besides the data generation and augmentation, our proposed approaches also include novel ML-driven policies for personalized solver strategies, with an emphasis on applications like graph convolutional networks for initial basis selection and reinforcement learning for advanced presolving and cut selection. Additionally, we detail the incorporation of state-of-the-art parameter tuning algorithms which markedly elevate solver performance. Compared with traditional solvers such as Cplex and SCIP, our ML-augmented OptVerse AI Solver demonstrates superior speed and precision across both established benchmarks and real-world scenarios, reinforcing the practical imperative and effectiveness of machine learning techniques in mathematical programming solvers.

In this paper, we propose an efficient sieving based secant method to address the computational challenges of solving sparse optimization problems with least-squares constraints. A level-set method has been introduced in [X. Li, D.F. Sun, and K.-C. Toh, SIAM J. Optim., 28 (2018), pp. 1842--1866] that solves these problems by using the bisection method to find a root of a univariate nonsmooth equation $\varphi(\lambda) = \varrho$ for some $\varrho > 0$, where $\varphi(\cdot)$ is the value function computed by a solution of the corresponding regularized least-squares optimization problem. When the objective function in the constrained problem is a polyhedral gauge function, we prove that (i) for any positive integer $k$, $\varphi(\cdot)$ is piecewise $C^k$ in an open interval containing the solution $\lambda^*$ to the equation $\varphi(\lambda) = \varrho$; (ii) the Clarke Jacobian of $\varphi(\cdot)$ is always positive. These results allow us to establish the essential ingredients of the fast convergence rates of the secant method. Moreover, an adaptive sieving technique is incorporated into the secant method to effectively reduce the dimension of the level-set subproblems for computing the value of $\varphi(\cdot)$. The high efficiency of the proposed algorithm is demonstrated by extensive numerical results.

In this paper, we present a method to certify the approximation quality of a low rank tensor to a given third order symmetric tensor. Under mild assumptions, best low rank approximation is attained if a control parameter is zero or quantified quasi-optimal low rank approximation is obtained if the control parameter is positive.This is based on a primal-dual method for computing a low rank approximation for a given tensor. The certification is derived from the global optimality of the primal and dual problems, and is characterized by easily checkable relations between the primal and the dual solutions together with another rank condition. The theory is verified theoretically for orthogonally decomposable tensors as well as numerically through examples in the general case., Comment: 46pages

In this paper, we propose an adaptive sieving (AS) strategy for solving general sparse machine learning models by effectively exploring the intrinsic sparsity of the solutions, wherein only a sequence of reduced problems with much smaller sizes need to be solved. We further apply the proposed AS strategy to generate solution paths for large-scale sparse optimization problems efficiently. We establish the theoretical guarantees for the proposed AS strategy including its finite termination property. Extensive numerical experiments are presented in this paper to demonstrate the effectiveness and flexibility of the AS strategy to solve large-scale machine learning models.

The exclusive lasso (also known as elitist lasso) regularizer has become popular recently due to its superior performance on intra-group feature selection. Its complex nature poses difficulties for the computation of high-dimensional machine learning models involving such a regularizer. In this paper, we propose a highly efficient dual Newton method based proximal point algorithm (PPDNA) for solving large-scale exclusive lasso models. As important ingredients, we systematically study the proximal mapping of the weighted exclusive lasso regularizer and the corresponding generalized Jacobian. These results also make popular first-order algorithms for solving exclusive lasso models more practical. Extensive numerical results are presented to demonstrate the superior performance of the PPDNA against other popular numerical algorithms for solving the exclusive lasso problems., Comment: arXiv admin note: substantial text overlap with arXiv:2009.08719

To ensure the system stability of the $\bf{\mathcal{H}_{2}}$-guaranteed cost optimal decentralized control problem (ODC), an approximate semidefinite programming (SDP) problem is formulated based on the sparsity of the gain matrix of the decentralized controller. To reduce data storage and improve computational efficiency, the SDP problem is vectorized into a conic programming (CP) problem using the Kronecker product. Then, a proximal alternating direction method of multipliers (PADMM) is proposed to solve the dual of the resulted CP. By linking the (generalized) PADMM with the (relaxed) proximal point algorithm, we are able to accelerate the proposed PADMM via the Halpern fixed-point iterative scheme. This results in a fast convergence rate for the Karush-Kuhn-Tucker (KKT) residual along the sequence generated by the accelerated algorithm. Numerical experiments further demonstrate that the accelerated PADMM outperforms both the well-known CVXOPT and SCS algorithms for solving the large-scale CP problems arising from $\bf{\mathcal{H}_{2}}$-guaranteed cost ODC problems.

In this paper, we propose a randomly projected convex clustering model for clustering a collection of $n$ high dimensional data points in $\mathbb{R}^d$ with $K$ hidden clusters. Compared to the convex clustering model for clustering original data with dimension $d$, we prove that, under some mild conditions, the perfect recovery of the cluster membership assignments of the convex clustering model, if exists, can be preserved by the randomly projected convex clustering model with embedding dimension $m = O(\epsilon^{-2}\log(n))$, where $0 < \epsilon < 1$ is some given parameter. We further prove that the embedding dimension can be improved to be $O(\epsilon^{-2}\log(K))$, which is independent of the number of data points. Extensive numerical experiment results will be presented in this paper to demonstrate the robustness and superior performance of the randomly projected convex clustering model. The numerical results presented in this paper also demonstrate that the randomly projected convex clustering model can outperform the randomly projected K-means model in practice.

In this paper, we study variational inequality problems (VIPs) with involved mappings and feasible sets characterized by polynomial functions (namely, polynomial VIPs). We propose a numerical algorithm for computing solutions to polynomial VIPs based on Lagrange multiplier expressions and the Moment-SOS hierarchy of semidefinite relaxations. We also extend our approach to finding more or even all solutions to polynomial VIPs. We show that the method proposed in this paper can find solutions or detect the nonexistence of solutions within finitely many steps, under some general assumptions. In addition, we show that if the VIP is given by generic polynomials, then it has finitely many Karush-Kuhn-Tucker points, and our method can solve it within finitely many steps. Numerical experiments are conducted to illustrate the efficiency of the proposed methods.

This paper proposes a squared smoothing Newton method via the Huber smoothing function for solving semidefinite programming problems (SDPs). We first study the fundamental properties of the matrix-valued mapping defined upon the Huber function. Using these results and existing ones in the literature, we then conduct rigorous convergence analysis and establish convergence properties for the proposed algorithm. In particular, we show that the proposed method is well-defined and admits global convergence. Moreover, under suitable regularity conditions, i.e., the primal and dual constraint nondegenerate conditions, the proposed method is shown to have a superlinear convergence rate. To evaluate the practical performance of the algorithm, we conduct extensive numerical experiments for solving various classes of SDPs. Comparison with the state-of-the-art SDP solvers demonstrates that our method is also efficient for computing accurate solutions of SDPs., Comment: 49 pages

As a novel distributed learning paradigm, federated learning (FL) faces serious challenges in dealing with massive clients with heterogeneous data distribution and computation and communication resources. Various client-variance-reduction schemes and client sampling strategies have been respectively introduced to improve the robustness of FL. Among others, primal-dual algorithms such as the alternating direction of method multipliers (ADMM) have been found being resilient to data distribution and outperform most of the primal-only FL algorithms. However, the reason behind remains a mystery still. In this paper, we firstly reveal the fact that the federated ADMM is essentially a client-variance-reduced algorithm. While this explains the inherent robustness of federated ADMM, the vanilla version of it lacks the ability to be adaptive to the degree of client heterogeneity. Besides, the global model at the server under client sampling is biased which slows down the practical convergence. To go beyond ADMM, we propose a novel primal-dual FL algorithm, termed FedVRA, that allows one to adaptively control the variance-reduction level and biasness of the global model. In addition, FedVRA unifies several representative FL algorithms in the sense that they are either special instances of FedVRA or are close to it. Extensions of FedVRA to semi/un-supervised learning are also presented. Experiments based on (semi-)supervised image classification tasks demonstrate superiority of FedVRA over the existing schemes in learning scenarios with massive heterogeneous clients and client sampling.

In this paper, we propose and analyze an efficient Halpern-Peaceman-Rachford (HPR) algorithm for solving the Wasserstein barycenter problem (WBP) with fixed supports. While the Peaceman-Rachford (PR) splitting method itself may not be convergent for solving the WBP, the HPR algorithm can achieve an $O(1/\varepsilon)$ non-ergodic iteration complexity with respect to the Karush-Kuhn-Tucker (KKT) residual. More interestingly, we propose an efficient procedure with linear time computational complexity to solve the linear systems involved in the subproblems of the HPR algorithm. As a consequence, the HPR algorithm enjoys an $O({\rm Dim(P)}/\varepsilon)$ non-ergodic computational complexity in terms of flops for obtaining an $\varepsilon$-optimal solution measured by the KKT residual for the WBP, where ${\rm Dim(P)}$ is the dimension of the variable of the WBP. This is better than the best-known complexity bound for the WBP. Moreover, the extensive numerical results on both the synthetic and real data sets demonstrate the superior performance of the HPR algorithm for solving the large-scale WBP.

In this paper, we propose a novel adaptive sieving (AS) technique and an enhanced AS (EAS) technique, which are solver independent and could accelerate optimization algorithms for solving large scale convex optimization problems with intrinsic structured sparsity. We establish the finite convergence property of the AS technique and the EAS technique with inexact solutions of the reduced subproblems. As an important application, we apply the AS technique and the EAS technique on the convex clustering model, which could accelerate the state-of-the-art algorithm SSNAL by more than 7 times and the algorithm ADMM by more than 14 times.

Estimation of the precision matrix (or inverse covariance matrix) is of great importance in statistical data analysis and machine learning. However, as the number of parameters scales quadratically with the dimension $p$, computation becomes very challenging when $p$ is large. In this paper, we propose an adaptive sieving reduction algorithm to generate a solution path for the estimation of precision matrices under the $\ell_1$ penalized D-trace loss, with each subproblem being solved by a second-order algorithm. In each iteration of our algorithm, we are able to greatly reduce the number of variables in the {problem} based on the Karush-Kuhn-Tucker (KKT) conditions and the sparse structure of the estimated precision matrix in the previous iteration. As a result, our algorithm is capable of handling datasets with very high dimensions that may go beyond the capacity of the existing methods. Moreover, for the sub-problem in each iteration, other than solving the primal problem directly, we develop a semismooth Newton augmented Lagrangian algorithm with global linear convergence rate on the dual problem to improve the efficiency. Theoretical properties of our proposed algorithm have been established. In particular, we show that the convergence rate of our algorithm is asymptotically superlinear. The high efficiency and promising performance of our algorithm are illustrated via extensive simulation studies and real data applications, with comparison to several state-of-the-art solvers.

In this paper, we aim to solve high dimensional convex quadratic programming (QP) problems with a large number of quadratic terms, linear equality and inequality constraints. In order to solve the targeted {\bf QP} problems to a desired accuracy efficiently, we develop a two-phase {\bf P}roximal {\bf A}ugmented {\bf L}agrangian method {(QPPAL)}, with Phase I to generate a reasonably good initial point to warm start Phase II to obtain an accurate solution efficiently. More specifically, in Phase I, based on the recently developed symmetric Gauss-Seidel (sGS) decomposition technique, we design a novel sGS based semi-proximal augmented Lagrangian method for the purpose of finding a solution of low to medium accuracy. Then, in Phase II, a proximal augmented Lagrangian algorithm is proposed to obtain a more accurate solution efficiently. Extensive numerical results evaluating the performance of {QPPAL} against {existing state-of-the-art solvers Gurobi, OSQP and QPALM} are presented to demonstrate the high efficiency and robustness of our proposed algorithm for solving various classes of large-scale convex QP problems. {The MATLAB implementation of the software package QPPAL is available at: \url{https://blog.nus.edu.sg/mattohkc/softwares/qppal/}., Comment: 28 pages, 4 figures

Shape-constrained convex regression problem deals with fitting a convex function to the observed data, where additional constraints are imposed, such as component-wise monotonicity and uniform Lipschitz continuity. This paper provides a unified framework for computing the least squares estimator of a multivariate shape-constrained convex regression function in $\mathbb{R}^d$. We prove that the least squares estimator is computable via solving an essentially constrained convex quadratic programming (QP) problem with $(d+1)n$ variables, $n(n-1)$ linear inequality constraints and $n$ possibly non-polyhedral inequality constraints, where $n$ is the number of data points. To efficiently solve the generally very large-scale convex QP, we design a proximal augmented Lagrangian method (proxALM) whose subproblems are solved by the semismooth Newton method (SSN). To further accelerate the computation when $n$ is huge, we design a practical implementation of the constraint generation method such that each reduced problem is efficiently solved by our proposed proxALM. Comprehensive numerical experiments, including those in the pricing of basket options and estimation of production functions in economics, demonstrate that our proposed proxALM outperforms the state-of-the-art algorithms, and the proposed acceleration technique further shortens the computation time by a large margin., Comment: arXiv admin note: substantial text overlap with arXiv:2002.11410

We consider the problem of learning a graph under the Laplacian constraint with a non-convex penalty: minimax concave penalty (MCP). For solving the MCP penalized graphical model, we design an inexact proximal difference-of-convex algorithm (DCA) and prove its convergence to critical points. We note that each subproblem of the proximal DCA enjoys the nice property that the objective function in its dual problem is continuously differentiable with a semismooth gradient. Therefore, we apply an efficient semismooth Newton method to subproblems of the proximal DCA. Numerical experiments on various synthetic and real data sets demonstrate the effectiveness of the non-convex penalty MCP in promoting sparsity. Compared with the existing state-of-the-art method, our method is demonstrated to be more efficient and reliable for learning graph Laplacian with MCP., Comment: 32 pages

In this paper, we adopt the augmented Lagrangian method (ALM) to solve convex quadratic second-order cone programming problems (SOCPs). Fruitful results on the efficiency of the ALM have been established in the literature. Recently, it has been shown in [Cui, Sun, and Toh, {\em Math. Program.}, 178 (2019), pp. 381--415] that if the quadratic growth condition holds at an optimal solution for the dual problem, then the KKT residual converges to zero R-superlinearly when the ALM is applied to the primal problem. Moreover, Cui, Ding, and Zhao [{\em SIAM J. Optim.}, 27 (2017), pp. 2332-2355] provided sufficient conditions for the quadratic growth condition to hold under the metric subregularity and bounded linear regularity conditions for solving composite matrix optimization problems involving spectral functions. Here, we adopt these recent ideas to analyze the convergence properties of the ALM when applied to SOCPs. To the best of our knowledge, no similar work has been done for SOCPs so far. In our paper, we first provide sufficient conditions to ensure the quadratic growth condition for SOCPs. With these elegant theoretical guarantees, we then design an SOCP solver and apply it to solve various classes of SOCPs, such as minimal enclosing ball problems, classical trust-region subproblems, square-root Lasso problems, and DIMACS Challenge problems. Numerical results show that the proposed ALM based solver is efficient and robust compared to the existing highly developed solvers, such as Mosek and SDPT3., Comment: 25 pages, 0 figure

The doubly nonnegative (DNN) cone, being the set of all positive semidefinite matrices whose elements are nonnegative, is a popular approximation of the computationally intractable completely positive cone. The major difficulty for implementing a Newton-type method to compute the projection of a given large scale matrix onto the DNN cone lies in the possible failure of the constraint nondegeneracy, a generalization of the linear independence constraint qualification for nonlinear programming. Such a failure results in the singularity of the Jacobian of the nonsmooth equation representing the Karush-Kuhn-Tucker optimality condition that prevents the semismooth Newton-CG method from solving it with a desirable convergence rate. In this paper, we overcome the aforementioned difficulty by solving a sequence of better conditioned nonsmooth equations generated by the augmented Lagrangian method (ALM) instead of solving one above mentioned singular equation. By leveraging on the metric subregularity of the normal cone associated with the positive semidefinite cone, we derive sufficient conditions to ensure the dual quadratic growth condition of the underlying problem, which further leads to the asymptotically superlinear convergence of the proposed ALM. Numerical results on difficult randomly generated instances and from the semidefinite programming library are presented to demonstrate the efficiency of the algorithm for computing the DNN projection to a very high accuracy., Comment: 28 pages, 0 figure

The exclusive lasso (also known as elitist lasso) regularization has become popular recently due to its superior performance on structured sparsity. Its complex nature poses difficulties for the computation of high-dimensional machine learning models involving such a regularizer. In this paper, we propose an adaptive sieving (AS) strategy for generating solution paths of machine learning models with the exclusive lasso regularizer, wherein a sequence of reduced problems with much smaller sizes need to be solved. In order to solve these reduced problems, we propose a highly efficient dual Newton method based proximal point algorithm (PPDNA). As important ingredients, we systematically study the proximal mapping of the weighted exclusive lasso regularizer and the corresponding generalized Jacobian. These results also make popular first-order algorithms for solving exclusive lasso models practical. Various numerical experiments for the exclusive lasso models have demonstrated the effectiveness of the AS strategy for generating solution paths and the superior performance of the PPDNA.

Estimation of Gaussian graphical models is important in natural science when modeling the statistical relationships between variables in the form of a graph. The sparsity and clustering structure of the concentration matrix is enforced to reduce model complexity and describe inherent regularities. We propose a model to estimate the sparse Gaussian graphical models with hidden clustering structure, which also allows additional linear constraints to be imposed on the concentration matrix. We design an efficient two-phase algorithm for solving the proposed model. We develop a symmetric Gauss-Seidel based alternating direction method of the multipliers (sGS-ADMM) to generate an initial point to warm-start the second phase algorithm, which is a proximal augmented Lagrangian method (pALM), to get a solution with high accuracy. Numerical experiments on both synthetic data and real data demonstrate the good performance of our model, as well as the efficiency and robustness of our proposed algorithm.

Shape-constrained convex regression problem deals with fitting a convex function to the observed data, where additional constraints are imposed, such as component-wise monotonicity and uniform Lipschitz continuity. This paper provides a comprehensive mechanism for computing the least squares estimator of a multivariate shape-constrained convex regression function in $\mathbb{R}^d$. We prove that the least squares estimator is computable via solving a constrained convex quadratic programming (QP) problem with $(n+1)d$ variables and at least $n(n-1)$ linear inequality constraints, where $n$ is the number of data points. For solving the generally very large-scale convex QP, we design two efficient algorithms, one is the symmetric Gauss-Seidel based alternating direction method of multipliers ({\tt sGS-ADMM}), and the other is the proximal augmented Lagrangian method ({\tt pALM}) with the subproblems solved by the semismooth Newton method ({\tt SSN}). Comprehensive numerical experiments, including those in the pricing of basket options and estimation of production functions in economics, demonstrate that both of our proposed algorithms outperform the state-of-the-art algorithm. The {\tt pALM} is more efficient than the {\tt sGS-ADMM} but the latter has the advantage of being simpler to implement.

A majorized accelerated block coordinate descent (mABCD) method in Hilbert space is analyzed to solve a sparse PDE-constrained optimization problem via its dual. The finite element approximation method is investigated. The attractive $O(1/k^2)$ iteration complexity of {the mABCD} method for the dual objective function values can be achieved. Based on the convergence result, we prove the robustness with respect to the mesh size $h$ for the mABCD method by establishing that asymptotically the infinite dimensional ABCD method and finite dimensional discretizations have the same convergence property, and the number of iterations of mABCD method remains almost constant as the discretization is refined., Comment: arXiv admin note: substantial text overlap with arXiv:1709.00005, arXiv:1708.09094, arXiv:1709.09539

Undirected graphical models have been especially popular for learning the conditional independence structure among a large number of variables where the observations are drawn independently and identically from the same distribution. However, many modern statistical problems would involve categorical data or time-varying data, which might follow different but related underlying distributions. In order to learn a collection of related graphical models simultaneously, various joint graphical models inducing sparsity in graphs and similarity across graphs have been proposed. In this paper, we aim to propose an implementable proximal point dual Newton algorithm (PPDNA) for solving the group graphical Lasso model, which encourages a shared pattern of sparsity across graphs. Though the group graphical Lasso regularizer is non-polyhedral, the asymptotic superlinear convergence of our proposed method PPDNA can be obtained by leveraging on the local Lipschitz continuity of the Karush-Kuhn-Tucker solution mapping associated with the group graphical Lasso model. A variety of numerical experiments on real data sets illustrates that the PPDNA for solving the group graphical Lasso model can be highly efficient and robust., Comment: 24 pages

In this paper, we consider high-dimensional nonconvex square-root-loss regression problems and introduce a proximal majorization-minimization (PMM) algorithm for these problems. Our key idea for making the proposed PMM to be efficient is to develop a sparse semismooth Newton method to solve the corresponding subproblems. By using the Kurdyka-{\L}ojasiewicz property exhibited in the underlining problems, we prove that the PMM algorithm converges to a d-stationary point. We also analyze the oracle property of the initial subproblem used in our algorithm. Extensive numerical experiments are presented to demonstrate the high efficiency of the proposed PMM algorithm., Comment: 34 pages, 8 tables

Powerful interior-point methods (IPM) based commercial solvers, such as Gurobi and Mosek, have been hugely successful in solving large-scale linear programming (LP) problems. The high efficiency of these solvers depends critically on the sparsity of the problem data and advanced matrix factorization techniques. For a large scale LP problem with data matrix $A$ that is dense (possibly structured) or whose corresponding normal matrix $AA^T$ has a dense Cholesky factor (even with re-ordering), these solvers may require excessive computational cost and/or extremely heavy memory usage in each interior-point iteration. Unfortunately, the natural remedy, i.e., the use of iterative methods based IPM solvers, although can avoid the explicit computation of the coefficient matrix and its factorization, is not practically viable due to the inherent extreme ill-conditioning of the large scale normal equation arising in each interior-point iteration. To provide a better alternative choice for solving large scale LPs with dense data or requiring expensive factorization of its normal equation, we propose a semismooth Newton based inexact proximal augmented Lagrangian ({\sc Snipal}) method. Different from classical IPMs, in each iteration of {\sc Snipal}, iterative methods can efficiently be used to solve simpler yet better conditioned semismooth Newton linear systems. Moreover, {\sc Snipal} not only enjoys a fast asymptotic superlinear convergence but is also proven to enjoy a finite termination property. Numerical comparisons with Gurobi have demonstrated encouraging potential of {\sc Snipal} for handling large-scale LP problems where the constraint matrix $A$ has a dense representation or $AA^T$ has a dense factorization even with an appropriate re-ordering., Comment: Due to the limitation "The abstract field cannot be longer than 1,920 characters", the abstract appearing here is slightly shorter than that in the PDF file

Nowadays, analysing data from different classes or over a temporal grid has attracted a great deal of interest. As a result, various multiple graphical models for learning a collection of graphical models simultaneously have been derived by introducing sparsity in graphs and similarity across multiple graphs. This paper focuses on the fused multiple graphical Lasso model which encourages not only shared pattern of sparsity, but also shared values of edges across different graphs. For solving this model, we develop an efficient regularized proximal point algorithm, where the subproblem in each iteration of the algorithm is solved by a superlinearly convergent semismooth Newton method. To implement the semismooth Newton method, we derive an explicit expression for the generalized Jacobian of the proximal mapping of the fused multiple graphical Lasso regularizer. Unlike those widely used first order methods, our approach has heavily exploited the underlying second order information through the semismooth Newton method. This can not only accelerate the convergence of the algorithm, but also improve its robustness. The efficiency and robustness of our proposed algorithm are demonstrated by comparing with some state-of-the-art methods on both synthetic and real data sets. Supplementary materials for this article are available online.

The exclusive lasso (also known as elitist lasso) regularization has become popular recently due to its superior performance on group sparsity. Compared to the group lasso regularization which enforces the competition on variables among different groups, the exclusive lasso regularization also enforces the competition within each group. In this paper, we propose a highly efficient dual Newton based preconditioned proximal point algorithm (PPDNA) to solve machine learning models involving the exclusive lasso regularizer. As an important ingredient, we provide a rigorous proof for deriving the closed-form solution to the proximal mapping of the weighted exclusive lasso regularizer. In addition, we derive the corresponding HS-Jacobian to the proximal mapping and analyze its structure --- which plays an essential role in the efficient computation of the PPA subproblem via applying a semismooth Newton method on its dual. Various numerical experiments in this paper demonstrate the superior performance of the proposed PPDNA against other state-of-the-art numerical algorithms.

This paper aims to present a fairly accessible generalization of several symmetric Gauss-Seidel decomposition based multi-block proximal alternating direction methods of multipliers (ADMMs) for convex composite optimization problems. The proposed method unifies and refines many constructive techniques that were separately developed for the computational efficiency of multi-block ADMM-type algorithms. Specifically, the majorized augmented Lagrangian functions, the indefinite proximal terms, the inexact symmetric Gauss-Seidel decomposition theorem, the tolerance criteria of approximately solving the subproblems, and the large dual step-lengths, are all incorporated in one algorithmic framework, which we named as sGS-imiPADMM. From the popularity of convergent variants of multi-block ADMMs in recent years, especially for high-dimensional multi-block convex composite conic programming problems, the unification presented in this paper, as well as the corresponding convergence results, may have the great potential of facilitating the implementation of many multi-block ADMMs in various problem settings.

We propose a semi-proximal augmented Lagrangian based decomposition method for convex composite quadratic conic programming problems with primal block angular structures. Using our algorithmic framework, we are able to naturally derive several well known augmented Lagrangian based decomposition methods for stochastic programming such as the diagonal quadratic approximation method of Mulvey and Ruszczy\'{n}ski. Moreover, we are able to derive novel enhancements and generalizations of these well known methods. We also propose a semi-proximal symmetric Gauss-Seidel based alternating direction method of multipliers for solving the corresponding dual problem. Numerical results show that our algorithms can perform well even for very large instances of primal block angular convex QP problems. For example, one instance with more than $300,000$ linear constraints and $12,500,000$ nonnegative variables is solved in less than a minute whereas Gurobi took more than 3 hours, and another instance {\tt qp-gridgen1} with more than $331,000$ linear constraints and $986,000$ nonnegative variables is solved in about 5 minutes whereas Gurobi took more than 35 minutes., Comment: 32 pages

In this paper, we study the polynomial optimization problem of multi-forms over the intersection of the multi-spheres and the nonnegative orthants. This class of problems is NP-hard in general, and includes the problem of finding the best nonnegative rank-one approximation of a given tensor. A Positivstellensatz is given for this class of polynomial optimization problems, based on which a globally convergent hierarchy of doubly nonnegative (DNN) relaxations is proposed. A (zero-th order) DNN relaxation method is applied to solve these problems, resulting in linear matrix optimization problems under both the positive semidefinite and nonnegative conic constraints. A worst case approximation bound is given for this relaxation method. Then, the recent solver SDPNAL+ is adopted to solve this class of matrix optimization problems. Typically, the DNN relaxations are tight, and hence the best nonnegative rank-one approximation of a tensor can be revealed frequently. Extensive numerical experiments show that this approach is quite promising., Comment: 27 pages

Spectral operators of matrices proposed recently in [C. Ding, D.F. Sun, J. Sun, and K.C. Toh, Math. Program. {\bf 168}, 509--531 (2018)] are a class of matrix valued functions, which map matrices to matrices by applying a vector-to-vector function to all eigenvalues/singular values of the underlying matrices. Spectral operators play a crucial role in the study of various applications involving matrices such as matrix optimization problems (MOPs) {that include semidefinite programming as one of the most important example classes}. In this paper, we will study more fundamental first- and second-order properties of spectral operators, including the Lipschitz continuity, $\rho$-order B(ouligand)-differentiability ($0<\rho\le 1$), $\rho$-order G-semismoothness ($0<\rho\le 1$), and characterization of generalized Jacobians., Comment: 25 pages. arXiv admin note: substantial text overlap with arXiv:1401.2269

Clustering is a fundamental problem in unsupervised learning. Popular methods like K-means, may suffer from poor performance as they are prone to get stuck in its local minima. Recently, the sum-of-norms (SON) model (also known as the clustering path) has been proposed in Pelckmans et al. (2005), Lindsten et al. (2011) and Hocking et al. (2011). The perfect recovery properties of the convex clustering model with uniformly weighted all pairwise-differences regularization have been proved by Zhu et al. (2014) and Panahi et al. (2017). However, no theoretical guarantee has been established for the general weighted convex clustering model, where better empirical results have been observed. In the numerical optimization aspect, although algorithms like the alternating direction method of multipliers (ADMM) and the alternating minimization algorithm (AMA) have been proposed to solve the convex clustering model (Chi and Lange, 2015), it still remains very challenging to solve large-scale problems. In this paper, we establish sufficient conditions for the perfect recovery guarantee of the general weighted convex clustering model, which include and improve existing theoretical results as special cases. In addition, we develop a semismooth Newton based augmented Lagrangian method for solving large-scale convex clustering problems. Extensive numerical experiments on both simulated and real data demonstrate that our algorithm is highly efficient and robust for solving large-scale problems. Moreover, the numerical results also show the superior performance and scalability of our algorithm comparing to the existing first-order methods. In particular, our algorithm is able to solve a convex clustering problem with 200,000 points in $\mathbb{R}^3$ in about 6 minutes., Comment: arXiv admin note: substantial text overlap with arXiv:1802.07091

We consider the problem of computing a Wasserstein barycenter for a set of discrete probability distributions with finite supports, which finds many applications in areas such as statistics, machine learning and image processing. When the support points of the barycenter are pre-specified, this problem can be modeled as a linear programming (LP) problem whose size can be extremely large. To handle this large-scale LP, we analyse the structure of its dual problem, which is conceivably more tractable and can be reformulated as a well-structured convex problem with 3 kinds of block variables and a coupling linear equality constraint. We then adapt a symmetric Gauss-Seidel based alternating direction method of multipliers (sGS-ADMM) to solve the resulting dual problem and establish its global convergence and global linear convergence rate. As a critical component for efficient computation, we also show how all the subproblems involved can be solved exactly and efficiently. This makes our method suitable for computing a Wasserstein barycenter on a large-scale data set, without introducing an entropy regularization term as is commonly practiced. In addition, our sGS-ADMM can be used as a subroutine in an alternating minimization method to compute a barycenter when its support points are not pre-specified. Numerical results on synthetic data sets and image data sets demonstrate that our method is highly competitive for solving large-scale Wasserstein barycenter problems, in comparison to two existing representative methods and the commercial software Gurobi.