One-Rank Linear Transformations and Fejer-Type Methods: An Overview

Abstract

Subgradient methods are frequently used for optimization problems. However, subgradient techniques are characterized by slow convergence for minimizing ravine convex functions. To accelerate subgradient methods, special linear non-orthogonal transformations of the original space are used. This paper provides an overview of these transformations based on Shor’s original idea. Two one-rank linear transformations of Euclidean space are considered. These simple transformations form the basis of variable metric methods for convex minimization that have a natural geometric interpretation in the transformed space. Along with the space transformation, a search direction and a corresponding step size must be defined. Subgradient Fejer-type methods are analyzed to minimize convex functions, and Polyak step size is used for problems with a known optimal objective value. Convergence theorems are provided together with the results of numerical experiments. Directions for future research are discussed.

1. Introduction

In the last few years, algorithms for convex optimization have revolutionized algorithm design, both for discrete and continuous optimization problems [1,2,3,4,5]. Popular convex optimization methods are subgradient-type methods [6,7,8,9].

Subgradient-type methods are characterized by slow convergence when minimizing ravine convex functions. This is due to the fact that an antisubgradient forms with an angle close to $\pi /2$ directed toward a minimum point. The original way out of this situation was first realized and proposed by N.Z. Shor in [10] (see also [11]). It is based on the use of linear non-orthogonal space transformations (namely, the space dilation operator), which allow for changing angles between the subgradient and the direction of a minimum. In particular, this leads to subgradient-type methods in argument space transformed by a linear operator. That is why such methods can be considered as variable metric methods. For other approaches to minimize ravine convex functions, e.g., gradient averaging or the heavy ball method are used; see [1,4,8] and the references therein.

Historically, the first variable metric method was the one for minimization twice continuously differentiable functions proposed by Davidon in [12] and developed by Fletcher and Powell in [13]. It is based on the idea of the quadratic approximation of the function being minimized, and it actually imitates the Newton–Raphson method without explicitly calculating second derivatives [14,15]. As a result, the “quasi-Newton type method” name has been attached to it, and its interpretation as a variable metric method is rarely found in the literature. This is also due to writing form of the Davidon–Fletcher–Powell method, which uses a symmetric matrix correction, which does not contribute to interpretation of it as a method in transformed argument space. This circumstance also applies to other methods of the quasi-Newton type [8,16,17], as well as to a number of conjugate gradient methods [14]. Close to variable metric methods are mirror descent methods [18,19,20,21,22]. Note that the construction of fast and reliable mirror descent methods is an important problem.

Henceforth, variable metric methods will be considered to be those methods that use linear non-orthogonal space transformations, the convergence results of which are based on the study of the characteristic behavior of a function being minimized in the transformed space of arguments. The majority of variable metric methods are subgradient-type methods that use a space dilation operation both in the direction of the subgradient and in the direction of the difference of two consecutive subgradients (r-algorithms) [23,24,25,26]. The ellipsoid method also belongs to the first group, and the r-algorithms have proven themselves to be an effective tool for solving nonsmooth problems. Variable metric methods can also include two methods [27] which are essentially variants of r-algorithms with a variable coefficient of space dilation and classical Fejer adjustment of the step multiplier. A family of economical—in the sense of information processing—variable metric methods for convex programming problems allows one to obtain a fairly general scheme of the methods of simple-body centroids [28], which, along with space dilation, involves the use of other linear non-orthogonal space transformations.

Results of the convergence of the methods given above are based on the study of the behavior of the subgradient norm of a function being minimized in dilated space, the study of the behavior of the distance to the minimum point function in transformed space, and the use of monotonous reduction of the volume of the localization region of the extrema set. Obviously, there is a simple mechanism that allows one to extend the class of variable metric methods. This can be achieved both by increasing the number of linear non-orthogonal transformations and the number of investigated characteristics of a function being minimized in the transformed space of arguments. At the same time, preserving the gradient nature of the methods in a transformed space, it is possible to construct methods that are close to quasi-Newton methods in efficiency but that are applicable to a wider class of functions than smooth ones, including nonsmooth ones.

The practical effectiveness of r-algorithms for ravine problems is due to the fact that space dilation is aimed at improving the structure of the level surfaces of a function being minimized in the next space of arguments. However, in addition to space dilation, other linear non-orthogonal operators can be used for this. Two of these, which allow one to improve the structure of level surfaces in ravine functions, will be discussed in this paper. These transformations are quite simple, and for convex function minimization, it is possible to construct and substantiate variable metric methods that have a visual geometric interpretation in the transformed space of arguments. In addition to space transformation, the construction of such methods requires both the selection of movement direction from the current point and a certain method of adjusting the step multiplier in this direction. We will limit ourselves to the consideration of subgradient-type methods involving classical Fejer adjustment of the step multiplier in the direction of the antisubgradient in relation to convex function minimization with a known value of the function at the optimum point. This paper was prepared using References [29,30].

The content of this paper is presented as follows: Section 2 gives the formulation of the problem and briefly analyzes the Fejer-type methods. In Section 3, the one-rank ellipsoidal transformation of space is considered, and two simple methods with classical Fejer step control are given based on this transformation. In Section 4, the to-orthogonalizing space transformation and the family of orthogonal subgradient descent methods based on it with classical Fejer step control are considered. The practical effectiveness of the above methods is supported by numerical experiments. Section 3 and Section 4 also discuss possible schemes of a number of other algorithms, taking into account these transformations, as well as the constructive solution of a number of fundamental problems of Fejer-type methods for solving the convex programming problem.

2. Fejer-Type Methods

Let there be an unconditional minimization problem

$$minf\left(x\right),$$

(1)

where $f\left(x\right)$ is a convex function of a vector argument $x\in X$, $X=R^n$; $R^n$ is Euclidean n-dimensional space with scalar product $(x,y)$. Here, X is an initial argument space, and n is the size of the problem. Let the extrema set $X^{*}$ of the problem (1) be non-empty and the minimum value of the function $f\left(x\right)$ be known: $f^{*}=f\left(x^{*}\right)$, $x^{*}\in X^{*}$. Without loss of generality, for convenience of further considerations, let us assume that the set $X^{*}$ consists of the single point $x^{*}$.

This formulation of the problem has a number of applications. In particular, it is easy to reduce the problem of finding a feasible point of a consistent system of convex inequalities to the problem (1). If the system of convex inequalities is inconsistent, analysis of this problem simplifies obtaining a sufficient condition of inconsistency. In addition, application of the exact penalty functions method to solve smooth constrained convex minimization problems leads to nonsmooth convex minimization problems [8,31].

Since we will further discuss numerical methods for solving the problem (1), let us formulate what we will consider as its solution, or, in other words, a stopping criterion for these methods. Let $X_{\epsilon }^{*}=\left\{x:f\left(x\right)-f^{*}\le \epsilon \right\}$, where $\epsilon >0$. We will consider the point $x_k\in X_{\epsilon }^{*}$, for which $f\left(x_k\right)-f^{*}\le \epsilon$ is fulfilled, as a solution of the problem (1) with accuracy $\epsilon$ by functional ($\epsilon$-solution). The input parameter of the methods will be the accuracy ${\epsilon }_f$, with which the problem (1) is to be solved by a functional, and a stop is performed as soon as the point $x_k\in X_{{\epsilon }_f}^{*}$ is reached.

Fejer methods go back to [32,33] and belong to the subgradient descent type, which involves adjustment of a step multiplier in the direction of a normalized subgradient. In order to solve the problem (1), Fejer methods with Polyak step size are formulated using the following iterative process:

$$x_{k+1}=x_k-h_k\frac{\partial f\left(x_k\right)}{\parallel \partial f\left(x_k\right)\parallel },h_k=\gamma \frac{\left(f\left(x_i\right)-f^{*}\right)}{\parallel \partial f\left(x_k\right)\parallel },$$

(2)

where $\partial f\left(x_k\right)$ is a subgradient of $f\left(x\right)$ at the point $x_k$, and $\gamma$ is a scalar such that $0<\gamma <2$ [34]. This selection of the parameter $\gamma$ values guarantees a monotonous reduction of the distance to a minimum at each step of the process (2) and ensures the simplicity of proving the convergence of these processes. The $h_k$ step at $\gamma =1$ we will call the classical Fejer step, and the corresponding method we will call the classical Fejer method. For a modern presentation of the theory behind Fejer methods, see [5,35,36].

When solving the problem (1) with accuracy ${\epsilon }_f$ using a functional, the classical Fejer method is quite constructive and has the following form.

At the beginning of the process, we have $x_0\in R^n$, ${\epsilon }_f>0$.

Let $x_k\in R^n$, $f\left(x_k\right)$, and $\partial f\left(x_k\right)$ be obtained at the k-th step, these being the values of the function and subgradient calculated within it. If $f\left(x_k\right)-f^{*}\le {\epsilon }_f$, then $x_k$ is the sought point. Otherwise, we calculate another approximation

$$x_{k+1}=x_k-h_k\frac{\partial f\left(x_k\right)}{\parallel \partial f\left(x_k\right)\parallel },h_k=\frac{\left(f\left(x_i\right)-f^{*}\right)}{\parallel \partial f\left(x_k\right)\parallel }$$

(3)

and proceed to the $(k+1)$-th step.

The following theorem is valid.

Theorem 1

([9]). There exists a constant $C>0$ for which $\parallel \partial f\left(x_k\right)\parallel \le C$ for all steps of the method (3). Then

$$\underset{i=0,1,\dots ,k}{min}f\left(x_i\right)-f^{*}\le \frac{C\parallel x_0-x^{*}\parallel }{\sqrt{k+1}}$$

for any $k\ge 0$. Therefore, the method (3) allows one to find the ${\epsilon }_f$-solution of the problem (1) in no more than $K=\left[{\left(\frac{C\parallel x_0-x^{*}\parallel }{{\epsilon }_f}\right)}^2\right]+1$ steps, where $\left[a\right]$ is an integer part of the number a.

The proof of this theorem is based on the fact that at each step of the method (3), distance to the minimum point decreases with the accuracy of the shift square in the direction of the normalized subgradient as well as the fact that the step $h_k$ is limited by the size of the $({\epsilon }_f/C)$ value. We note that Theorem 1 condition can be strengthened if $\parallel \partial f\left(x\right)\parallel$ is bounded in a sphere of radius $\parallel x_0-x^{*}\parallel$ centered at the point $x^{*}$, $x_0\in R^n$.

Despite its theoretical appeal, the method in (3) has a bad reputation in practical applications even for smooth convex functions due to its slow convergence for ravine functions. In particular, when the level surfaces of a function $f\left(x\right)$ are not very dilated, it works well for even large-sized problems ($n\sim 100$). However, if the direction of the antisubgradient at the point is quite close to the orthogonal one (in the direction of the minimum, which is typical for ravine functions), then method (3) is almost hopeless even for small problems ($n\sim$ 10). For nonsmooth convex functions, where ravineness is the rule rather than an exception, the method (3) is unreliable for problems of even smaller sizes ($n\sim$ 2–5). In particular, minimizing two-variable function $f(x_1,x_2)=|x_1|+k|x_2|$ causes it to converge with the geometric progression rate with denominator $\sqrt{1-1/k^2}$, which is close to one for small k values ($k\sim$ 10–100).

To accelerate convergence, based on a piecewise linear approximation of a functional near the minimum, a more general method is proposed in [34]:

$$x_{k+1}=P_{Q_k}\left(x_k\right),Q_k=\left\{x:f\left(x_i\right)+(\partial f\left(x_i\right),x-x_i)\le f^{*},i\in I_k\right\}.$$

(4)

Here, $I_k$ is any subset of indices from $0,1,\dots ,k$ that necessarily contains k, and $P_{Q_k}$ is a projection operator on $Q_k$. The main problem of the method shown in (5) is to specify a constructive way of forming the set $I_k$ in order to limit the number of subgradients being stored. The most meaningful answer here the orthogonal descent method [37] can be considered. It is based on the idea of successive immersion of a non-obtuse cone containing $X^{*}$ into a rectangular cone, and uses no more than $(n+1)$ vectors. However, it does not always give good results in terms of the accuracy of solving the problem (5) by functional. This is evidenced by the results of a number of numerical experiments with both the orthogonal descent method and its modification, presented in [38,39]. Despite the small size of the problems ($n\sim 10$) and the fact that they are not very robust, the accuracy of their solution by functional is often ${\epsilon }_f\sim {10}^{-4}$–${10}^{-2}$.

Another Fejer-type method is proposed in [40]. It uses only two vectors, and to construct the next descent direction, it uses a linear combination of the subgradient direction and the previous step direction (like conjugate gradients) so that it forms a sharper angle with the direction to the minimum. For ravine functions, it provides a “smoother” movement along the ravine than the above methods using only two vectors, the behavior of which is “jumpy”. However, it is also inefficient for ravine functions, especially when the ravine is multidimensional.

In terms of the accuracy of solving the problem (1), to improve the situation for method (3), one can use linear transformation of space to align the structure of the level surfaces of the function in transformed space.

Let $Y=AX$ be the space of arguments transformed using the linear operator A. Here, A is a non-degenerate matrix of dimensions $n\times n$. Let B be the inverse operator of A: $B=A^{-1}$. Then, the classical Fejer method with the space transformation into X is as follows:

$$x_{k+1}=x_k-h_{k}B\frac{B^T\partial f\left(x_k\right)}{\parallel B^T\partial f\left(x_k\right)\parallel },h_k=\frac{\left(f\left(x_i\right)-f^{*}\right)}{\parallel B^T\partial f\left(x_k\right)\parallel },$$

(5)

which corresponds to the classical Fejer step in the direction of the antisubgradient in the transformed argument space Y

$$y_{k+1}=y_k-h_k\frac{\partial \phi \left(y_k\right)}{\parallel \partial \phi \left(y_k\right)\parallel },h_k=\frac{(\phi \left(y_k\right)-{\phi }^{*})}{\parallel \partial \phi \left(y_k\right)\parallel }.$$

(6)

Here, $y=Ax$ are images of points from X in the transformed space Y, and $\phi \left(y\right)$ is a convex function defined in the transformed space Y: $\phi \left(y\right)=f\left(A^{-1}y\right)=f\left(x\right)$. $\partial \phi \left(y_k\right)=B^T\partial f\left(x_k\right)$ is the subgradient at the point $y_k=A{x}_k$ for the function $\phi \left(y\right)$.

Let $Y_k=A_{k}X$ be the current argument space. Then, in terms of improving the structure of the level surfaces ${\phi }_{k+1}\left(y\right)$ in the next transformed space of arguments $Y_{k+1}=T_{k+1}{Y}_k$, it is natural to expect more efficient work from methods of type (5) than from method (3). Here, $T_{k+1}$ is a non-degenerate matrix of size $n\times n$ which specifies the transformation from $Y_k$ to $Y_{k+1}$.

The main relations that are needed in order to implement such a process are the recalculation of the matrix $B_{k+1}$ according the relation

$$B_{k+1}=B_k{T}_{k+1}^{-1}$$

(7)

and the recalculation of subgradients for ${\phi }_{k+1}\left(y\right)$ defined in $Y_{k+1}$ according to relation

$$\partial {\phi }_{k+1}\left(y\right)={\left(T_{k+1}^{-1}\right)}^T\partial {\phi }_k\left(y\right).$$

(8)

Further, we will often use various identities that follow from the chain of equalities

$$A_{k+1}=B_{k+1}^{-1}={\left(B_k{T}_{k+1}^{-1}\right)}^{-1}=T_{k+1}{B}_k^{-1}=T_{k+1}{A}_k.$$

Taking into account

$$\parallel B^T\partial f\left(x_k\right){\parallel }^2=\left(B^T\partial f\left(x_k\right),B^T\partial f\left(x_k\right)\right)=\left(\partial f\left(x_k\right),B{B}^T\partial f\left(x_k\right)\right),$$

process (5) can be written in the form

$$x_{k+1}=x_k-\frac{\left(f\left(x_k\right)-f^{*}\right)}{\left(\partial f\left(x_k\right),H\partial f\left(x_k\right)\right)}H\partial f\left(x_k\right),$$

(9)

where $H=B{B}^T$ is a positive definite symmetric matrix of dimension $n\times n$. Depending on which of the ratios are taken as a basis, there are two means of practical implementation of the methods. When referring to them, we will stick to the names that have been fixed for them for the subgradient-type methods with space dilation. The method based on (5) is usually called the method in B form (the B method), and the method based on (9) is called the method in the H form (the H method). The advantage of the H method is that for storing the symmetric matrix H, it requires almost twice as much RAM ($n\times (n+1)/2$) as the B method ($n\times n$ cells). However, the method in the B form allows for simple interpretation of the process in the transformed space Y. Therefore, it is advisable to build methods in the B form, and, as a consequence, to consider their analogues in the H form, mainly to save memory and calculations. Therefore, we will focus on methods based on ratio (5).

3. One-Rank Ellipsoidal Operator

To solve the problem (1) in [27], two subgradient methods with a classical Fejer step are considered in a transformed argument space. They are based on volume reduction of the localization region of the extrema set and preserve the monotony of estimation of the distance to the optimum in the transformed space of arguments. The space transformation used in them is similar to that used in r-algorithms [23,24,25], and is reduced to two consecutive dilations of the argument space in orthogonal directions. The space transformation operator has the form

$$T_2(\xi ,\eta )=R_{{\alpha }_1}\left(\frac{\xi +\eta }{\parallel \xi +\eta \parallel }\right)R_{{\alpha }_2}\left(\frac{\xi -\eta }{\parallel \xi -\eta \parallel }\right).$$

(10)

Here, $\xi$, $\eta$ are normalized vectors from $R^n$ such that ${(\xi ,\eta )}^2\ne 1$. ${\alpha }_1=1/{\beta }_1$, where ${\beta }_1=\sqrt{1-(\xi ,\eta )}$, and ${\alpha }_2=1/{\beta }_2$, where ${\beta }_2=\sqrt{1+(\xi ,\eta )}$. Further, we will call the operator (10) a two-rank ellipsoidal operator.

The methods in [27] use trivial information about hyperplanes that cut off the extrema set. In particular, the first method uses only two subgradients at consecutive points obtained according to the classical Fejer step. The second method, in addition to two consecutive subgradients, also uses a third vector of the aggregate type, which is a convex combination of previously calculated subgradients. Despite modesty in the information used, computational experiments using these methods described in [27] for a number of both smooth and nonsmooth test problems turned out to be not so bad. Therefore, more appropriate strategies for using calculated subgradients can lead to more efficient methods based on the transformation (10).

However, despite the simple geometric interpretation, the two-rank ellipsoidal transformation is inconvenient because, for methods in the B form, when moving from the k-th to the ($k+1$)-th step, the matrix $B_{k+1}$ calculation requires two one-rank corrections of a matrix of dimensions $n\times n$, resulting in $4n^2$ arithmetic operations. However, it turns out that there is a space transformation that allows for one-rank correction of the $B_{k+1}$ matrix and that basically does not differ from the two-rank ellipsoidal transformation. This is not surprising, since methods like (5) are defined by the matrix $H_k=B_k{B}_k^T$. Due to the ambiguity of the decomposition of matrix $H_k$, the same method in the H form may correspond to different methods in the B form.

Let $\xi$, $\eta \in R^n$ be vectors such that $\parallel \xi \parallel =1$, $\parallel \eta \parallel =1$, and their scalar product satisfies the condition ${(\xi ,\eta )}^2\ne 1$. A linear operator from $R^n$ to $R^n$ is called a one-rank ellipsoidal operator, which is written in the following matrix form:

$$T_1(\xi ,\eta )=I-\frac{1}{1-{(\xi ,\eta )}^2}\left(\left(1-\sqrt{1-{(\xi ,\eta )}^2}\right)\eta -(\xi ,\eta )\xi \right){\eta }^T.$$

(11)

Here, I is an unitary matrix of size $n\times n$. This operator has a number of interesting properties, which allows one to build and justify various variable metric methods on its basis. The following characteristics are valid for this matrix:

Lemma 1

([29]). For the operator $T_1(\xi ,\eta )$ given according to (11), if ${(\xi ,\eta )}^2\ne 1$, then there exists an inverse $T_1^{-1}(\xi ,\eta )$ and

$$T_1^{-1}(\xi ,\eta )=I+\frac{1}{\sqrt{1-{(\xi ,\eta )}^2}}\left(\left(1-\sqrt{1-{(\xi ,\eta )}^2}\right)\eta -(\xi ,\eta )\xi \right){\eta }^T.$$

(12)

In addition, for $T_1(\xi ,\eta )$ and $T_1^{-1}(\xi ,\eta )$, the following relations are fulfilled:

$$\begin{array}{ccc}\hfill T_1^T(\xi ,\eta )T_1(\xi ,\eta )& =& I+\frac{(\xi ,\eta )}{1-{(\xi ,\eta )}^2}\left(\xi {\eta }^T+\eta {\xi }^T\right),\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill T_1^{-1}(\xi ,\eta ){\left(T_1^{-1}(\xi ,\eta )\right)}^T& =& I+\frac{{(\xi ,\eta )}^2}{1-{(\xi ,\eta )}^2}\left(\xi {\xi }^T+\eta {\eta }^T\right)\hfill \\ & -& \frac{(\xi ,\eta )}{1-{(\xi ,\eta )}^2}\left(\xi {\eta }^T+\eta {\xi }^T\right).\hfill \end{array}$$

(14)

A proof for Lemma 1 is not presented due to its bulkiness. We only note that validity of (12)–(14) is easy to check by direct verification. An acceptable type of expressions (13) and (14) is caused by the fact that when summing similar terms with one-rank matrices $\xi {\xi }^T$, $\xi {\eta }^T$, $\eta {\xi }^T$, and $\eta {\eta }^T$, members with $\sqrt{1-{(\xi ,\eta )}^2}$ are reduced.

Lemma 2

([29]). Let $B_k$ be a non-degenerate matrix of dimensions $n\times n$, $p_1$, $p_2$—n-dimensional vectors such that ${\left(\frac{B_k^T{p}_1}{\parallel B_k^T{p}_1\parallel },\frac{B_k^T{p}_2}{\parallel B_k^T{p}_2\parallel }\right)}^2<1$. Let

$$B_{k+1}=B_k{T}_1^{-1}(\xi ,\eta ),\mathrm{where}\xi =\frac{B_k^T{p}_1}{\parallel B_k^T{p}_1\parallel }\mathrm{and}\eta =\frac{B_k^T{p}_2}{\parallel B_k^T{p}_2\parallel }.$$

(15)

Then, the matrix $B_{k+1}$ is non-degenerate, and the following properties are fulfilled for it:

$$(\mathbf{a}) det\left(B_{k+1}\right)=det\left(B_k\right)\sqrt{1-{(\xi ,\eta )}^2}\hspace{1em}and\hspace{1em}(\mathbf{b}) \left(B_{k+1}^T{p}_1,B_{k+1}^T{p}_2\right)=0$$

To prove Lemma 2, expression (14) from Lemma 1 and some technical features are used.

One-rank ellipsoidal operator (11) has a clear geometric interpretation (see Figure 1). It transforms ellipsoid $Ell(x_0,r)$ into a ball. The minimum volume ellipsoid $Ell(x_0,r)$ is centered at $x_0$ and contains the convex set resulting from the intersection of the ball $S(x_0,r):=\{x\in {\mathbb{R}}^n:\parallel x-x_0\parallel \le r\}$ with two half-spaces

$$\begin{array}{c}P(x_0,\xi ):=\{x\in {\mathbb{R}}^n:{(x-x_0)}^T\xi \le 0\},\\ P(x_0,\eta ):=\{x\in {\mathbb{R}}^n:{(x-x_0)}^T\eta \le 0\},\end{array}$$

where $-1<{\xi }^T\eta <0$, $\parallel \xi \parallel =1$, $\parallel \eta \parallel =1$. Thus, the vector $\xi$ remains unchangeable under transformation (11), while the vector $\eta$ is transformed so that $(\xi ,\eta )=0$ (see Figure 1b).

One-rank ellipsoidal space transformation (11) is strongly related to two-rank ellipsoidal transformation (10). Specifically, for both of them, the same decrease in inverse matrix determinant is characteristic, which is equivalent to a volume decrease of the localization area of the extrema set under certain conditions. Both the first and the second orthogonalize vectors when transitioning to the transformed space of arguments. Further, we will see that both of these transformations lead to the same method in the H form. However, unlike the two-rank ellipsoidal transformation, the transformation (11) requires half as many arithmetic operations when recalculating the matrix $B_k$, i.e., exactly as much as the space dilation in some direction. Therefore, its use leads to computationally more economical methods in the B form than the use of the two-rank ellipsoidal transformation (10).

Based on the volume reduction of the localization region of the extrema set, the one-rank ellipsoidal operator (11) simplifies, justifying methods for solving the problem (1) to improve the structure of level surfaces of the function in the next transformed space. Let $g_k=B_k^T\partial f\left(x_k\right)$ and $g_{k+1}=B_k^T\partial f\left(x_{k+1}\right)$ be subgradients of ${\phi }_k\left(y\right)$ in $Y_k=A_{k}X$ at the points $y_k=A_k{x}_k$ and $y_{k+1}=A_k{x}_{k+1}$. Here, $y_{k+1}$ is obtained according to the classic Fejer step in the transformed argument space $Y_k$. Then, if $(g_k,g_{k+1})\ge 0$, then the space transformation is not necessary. Let $(g_k,g_{k+1})<0$. Then, the space transformation $Y_{k+1}=T_1\left(\frac{g_k}{\parallel g_k\parallel },\frac{g_{k+1}}{\parallel g_{k+1}\parallel }\right)Y_k$ or $Y_{k+1}=T_1\left(\frac{g_{k+1}}{\parallel g_{k+1}\parallel },\frac{g_k}{\parallel g_k\parallel }\right)Y_k$ according to Lemma 2 allows one to orthogonalize the images of subgradients $g_k$ and $g_{k+1}$ in the transformed space of arguments $Y_{k+1}$ while providing better level surfaces for the function ${\phi }_{k+1}\left(y\right)$ in $Y_{k+1}$. These simple considerations lead to the following method for solving the problem (1).

3.1. Fejer-Type Method with One-Rank Ellipsoidal Space Transformation Using Two Successive Subgradients

Before starting calculations, we have $x_0\in R^n$; $B_0=I$ is a unit matrix of dimensions $n\times n$; ${\epsilon }_f$ is the accuracy of the functional with which the problem (1) should be solved; $f\left(x_0\right)$ and $g\left(x_0\right)=\partial f\left(x_0\right)$ are the function and subgradient values calculated in $x_0$. Then, if $f\left(x_0\right)-f^{*}\le {\epsilon }_f$, then $x_0$ is the sought point and STOP.

Let $x_k\in R^n$ be obtained at the iteration k, $f\left(x_k\right)$, $g\left(x_k\right)=\partial f\left(x_k\right)$, and $B_k$ be an $n\times n$-dimensional matrix.

Step 1. Calculate the next approximation

$$x_{k+1}=x_k-\frac{f\left(x_k\right)-f^{*}}{\parallel B_k^{T}g\left(x_k\right)\parallel }{B}_k\frac{B_k^{T}g\left(x_k\right)}{\parallel B_k^{T}g\left(x_k\right)\parallel }.$$

(16)

Step 2. Calculate $f\left(x_{k+1}\right)$ and $g\left(x_{k+1}\right)=\partial f\left(x_{k+1}\right)$. If $f\left(x_{k+1}\right)-f^{*}\le {\epsilon }_f$, then $x_{k+1}$ is the sought point and STOP. Otherwise, proceed to the Step 3.

Step 3. Set

$${\xi }_1=\frac{B_k^{T}g\left(x_k\right)}{\parallel B_k^{T}g\left(x_k\right)\parallel },{\xi }_2=\frac{B_k^{T}g\left(x_{k+1}\right)}{\parallel B_k^{T}g\left(x_{k+1}\right)\parallel }.$$

(17)

Step 4. If $({\xi }_1,{\xi }_2)\ge 0$, consider $B_{k+1}=B_k$ and go to the Step 5. Otherwise, calculate

$$B_{k+1}=B_k{T}_1^{-1}({\xi }_1,{\xi }_2)\mathrm{or}{B}_{k+1}=B_k{T}_1^{-1}({\xi }_2,{\xi }_1).$$

(18)

Step 5. Go to the iteration $(k+1)$ with $x_{k+1}$, $B_{k+1}$, $g\left(x_{k+1}\right)$, $f\left(x_{k+1}\right)$.

The following theorem holds.

Theorem 2

([29]). The sequence ${\left\{x_{k+1}\right\}}_{k=0}^{\infty }$, which is generated by the method (16)–(18), satisfies the inequalities

$$\parallel A_{k+1}(x_{k+1}-x^{*}){\parallel }^2\le {\parallel A_k(x_k-x^{*})\parallel }^2-\frac{{\left(f\left(x_k\right)-f^{*}\right)}^2}{\parallel B_k^{T}g\left(x_k\right){\parallel }^2}.$$

(19)

Here, $A_k=B_k^{-1}$, $A_{k+1}=B_{k+1}^{-1}$, $k=0,1,2,\dots$.

The proof of Theorem 2 is rather bulky and it uses statements from Lemma 1. For the method (16)–(18), Theorem 2 provides the possibility of substantiating its marginal convergence to the solution of problem (1) by reducing the ellipsoid volume that localizes $x^{*}$. At the same time, the volume reduction is better the larger the classical Fejer step in transformed space and the more obtuse the angle between successive subgradients. However, both of these circumstances depend on the specific properties of $f\left(x\right)$, and it is impossible to obtain any estimates for the general convexity of $f\left(x\right)$. Therefore, let us focus on this variant of proving the convergence of the method (16)–(18).

Theorem 3

([29]). Let the conditions $\parallel B_k\parallel \le c_1$ and $\parallel \partial f\left(x_k\right)\parallel \le c_2$ hold on each step the method (16)–(18). Here, $\parallel B_k\parallel$ is the Euclidean norm of the matrix $B_k$, i.e., $\parallel B_k\parallel ={\left({\displaystyle \sum _{i=1}^n}{\displaystyle \sum _{j=1}^n}{\left|b_{ij}^k\right|}^2\right)}^{1/2}$. Then, the method (16)–(18) solves the problem (1) with an accuracy ${\epsilon }_f$ by functional in no more than K steps, where $K=\left[{\left({\displaystyle \frac{c_1{c}_2\parallel x_0-x^{*}\parallel }{{\epsilon }_f}}\right)}^2\right]+1$. Here, $\left[a\right]$ is an integer part of the number a.

The proof of Theorem 3 is based on Theorem 1 and an assumption trick about the method (16)–(18).

Therefore, if the conditions of Theorem 3 are fulfilled, then the method (16)–(18) from any initial approximation converges to the ${\epsilon }_f$-solution of the problem (1) for a finite number of iterations. However, the estimate of the maximum number of iterations required to achieve ${\epsilon }_f$ accuracy is very rough. In real work using the method, it should be expected that for ravine functions, the stopping criterion will work with the condition that

$$f\left(x_k\right)-f^{*}\le \parallel B_k^{T}g\left(x_k\right)\parallel \sqrt{\parallel x_0-x^{*}{\parallel }^2-\sum _{i=0}^{k-1}{\displaystyle \frac{{\left(f\left(x_i\right)-f^{*}\right)}^2}{\parallel B_i^{T}g\left(x_i\right){\parallel }^2}}},$$

(20)

for the value of $\parallel B_k^{T}g\left(x_k\right)\parallel$ will be rather small for large k due to the fact that $det\left(B_k\right)$ goes to zero as $k\to \infty$. Therefore, it is advisable to strengthen the space transformation in order to guarantee a stronger convergence to zero of $det\left(B_k\right)$ as $k\to \infty$. In particular, a variant of the method as in [27] is suitable for this, which uses an aggregate vector, a convex combination of previously calculated subgradients, to enhance the transformation. In this case, the method will become a little more complicated, and it, relative to the method (16)–(18), will be characterized by a stronger alignment of the level surfaces of the function in transformed space, especially in the case of multidimensional ravines. Let us focus on computational schemes of these two methods, which are convenient for implementation, although it should be noted that the transformation (11) allows for justifying a whole family of Fejer-type methods. We will discuss some of them in detail later.

Practical implementation of the methods. Despite the ease of proof, the method (16)–(18) is inconvenient to implement on a computer. In particular, when the space transformation operation is implemented, it requires calculation of the subgradient image $\partial f\left(x_{k+1}\right)$ as in the space $Y_k-B_k^T\partial f\left(x_{k+1}\right)$ as well as in the space $Y_{k+1}-B_{k+1}^T\partial f\left(x_{k+1}\right)$. Calculation of each of them requires $n^2$ arithmetic operations of multiplication and the same number of additions, which results in a total of $4n^2$ arithmetic operations. This can be avoided if the subgradient when moving to the next transformed argument space is recalculated according to the ratio (8). Moreover, in the methods considered the antisubgradient in the transformed space of arguments as the direction of movement. The following lemma allows this to be achieved for the transformation (11).

Lemma 3

([29]). Let $B_k$ be a non-degenerate matrix of size $n\times n$; $p_k$, $p_{k+1}$ are n-dimensional vectors such that ${\left({\displaystyle \frac{B_k^T{p}_k}{\parallel B_k^T{p}_k\parallel }},{\displaystyle \frac{B_k^T{p}_{k+1}}{\parallel B_k^T{p}_{k+1}\parallel }}\right)}^2<1$. Let $B_{k+1}=T_1^{-1}({\xi }_k,{\xi }_{k+1})$, where ${\xi }_k={\displaystyle \frac{B_k^T{p}_k}{\parallel B_k^T{p}_k\parallel }}$, ${\xi }_{k+1}={\displaystyle \frac{B_k^T{p}_{k+1}}{\parallel B_k^T{p}_{k+1}\parallel }}$. Then,

$$\begin{array}{ccc}\hfill {\displaystyle \frac{B_{k+1}^T{p}_{k+1}}{\parallel B_{k+1}^T{p}_{k+1}\parallel }}& =& {\xi }_{k+1},\hfill \end{array}$$

$$\begin{array}{ccc}\hfill {\displaystyle \frac{B_{k+1}^T{p}_k}{\parallel B_{k+1}^T{p}_k\parallel }}& =& {\displaystyle \frac{1}{\sqrt{1-{({\xi }_k,{\xi }_{k+1})}^2}}}\left({\xi }_k-({\xi }_k,{\xi }_{k+1}){\xi }_{k+1}\right),\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill \parallel B_{k+1}^T{p}_i\parallel & =& \sqrt{1-{({\xi }_k,{\xi }_{k+1})}^2}\parallel B_k^T{p}_i\parallel ,i=k,k+1.\hfill \end{array}$$

(22)

For the vectors $B_k^T{p}_k$ and $B_k^T{p}_{k+1}$ in transformed space $Y_k=A_{k}X$, the expressions (21) and (22) of Lemma 3 ensure recalculation of both normalized directions and their norms at the transition to the next transformed space of arguments obtained as a result of applying the operator $T_1\left(\frac{B_k^T{p}_k}{\parallel B_k^T{p}_k\parallel },\frac{B_k^T{p}_{k+1}}{|B_k^T{p}_{k+1}\parallel }\right)$. However, for the B form of the methods, it does not matter which of the transformations—$T_1\left(\frac{B_k^T{p}_k}{\parallel B_k^T{p}_k\parallel },\frac{B_k^T{p}_{k+1}}{\parallel B_k^T{p}_{k+1}\parallel }\right)$ or $T_1\left(\frac{B_k^T{p}_{k+1}}{\parallel B_k^T{p}_{k+1}\parallel },\frac{B_k^T{p}_k}{\parallel B_k^T{p}_k\parallel }\right)$—is chosen. The choice of the first one is more rational, since it ensures recalculation of subgradients, which allows one not to accumulate errors in normalized directions for making a step in the next transformed space of arguments.

For the implementation of the method, recalculation of the classic Fejer step multiplier in the direction of the normalized subgradient is also required during the transition to the next transformed space of arguments. This does not lead to any problems, since their recalculation is associated with changing the norms of subgradients and, therefore,

$$h_{k+1}=h_k\frac{\parallel B_k^T{g}_k\parallel }{\parallel B_{k+1}^T{g}_k\parallel }.$$

Here, $h_{k+1}$ and $h_k$ are classical Fejer steps in the direction of the normalized subgradient in the spaces $Y_{k+1}=A_{k+1}X$ and $Y_k=A_{k}X$.

Therefore, Lemma 3 allows one to avoid unnecessary calculations for the B form of the method (16)–(18), and for the problem (1), it takes on a more rational form in terms of computation.

Before starting calculations, we have ${\epsilon }_f>0$, $x_0\in R^n$. Then, if $f\left(x_0\right)-f^{*}\le {\epsilon }_f$, then $x_0$ is the sought point and STOP. Otherwise, we set $h_0=\frac{f\left(x_0\right)-f^{*}}{\parallel \partial f\left(x_0\right)\parallel }$, ${\xi }_0=\frac{\partial f\left(x_0\right)}{\parallel \partial f\left(x_0\right)\parallel }\in R^n$; $B_0=I$ is a unit matrix of dimension $n\times n$.

Let the $x_k\in R^n$, $h_k$, ${\xi }_k\in R^n$; $B_k$ matrix of dimension $n\times n$ be obtained at the iteration k.

Step 1. Calculate the next approximation

$$\begin{array}{cc}\hfill & x_{k+1}=x_k-h_k{B}_k{\xi }_k.\end{array}$$

(23)

Step 2. Calculate $f\left(x_{k+1}\right)$, $\partial f\left(x_{k+1}\right)$. If $f\left(x_{k+1}\right)-f^{*}\le {\epsilon }_f$, then $x_{k+1}$ is the sought point and STOP. Otherwise, set

$${\xi }_{k+1}={\displaystyle \frac{B_k^T\partial f\left(x_{k+1}\right)}{\parallel B_k^T\partial f\left(x_{k+1}\right)|}},h_{k+1}={\displaystyle \frac{f\left(x_{k+1}\right)-f^{*}}{\parallel B_k^T\partial f\left(x_{k+1}\right)\parallel }}.$$

(24)

Step 3. If $({\xi }_k,{\xi }_{k+1})\ge 0$, then $B_{k+1}=B_k$ and go to the Step 4. Otherwise, calculate

$$\eta =\left({\displaystyle \frac{1}{\sqrt{1-{({\xi }_k,{\xi }_{k+1})}^2}}}-1\right){\xi }_{k+1}-{\displaystyle \frac{({\xi }_k,{\xi }_{k+1})}{\sqrt{1-{({\xi }_k,{\xi }_{k+1})}^2}}}{\xi }_k,$$

(25)

$$B_{k+1}=B_k\left(I+\eta {\xi }_{k+1}^T\right),h_{k+1}={\displaystyle \frac{h_{k+1}}{\sqrt{1-{({\xi }_k,{\xi }_{k+1})}^2}}}.$$

(26)

Step 4. Go to the iteration $(k+1)$ with $x_{k+1}$, $B_{k+1}$, ${\xi }_{k+1}$, $h_{k+1}$.

Therefore, the method (23)–(26) is identical to the method (16)–(18), and Theorems 2 and 3 hold for it. However, it is computationally more rational than the method (23)–(26), and further, we will choose it for numerical experiments.

3.2. A Fejer-Type Method with One-Rank Ellipsoidal Space Transformation Using Two Consecutive Subgradients and an Aggregate-Type Vector

Identically to [27], it is easy to construct a method that uses two consecutive subgradients and some aggregate vector that allows us to strengthen the space transformation and ensure faster convergence of the method in the sense of (20). Using Lemmas 2 and 3 given above, as well as Lemma 3 from [27], we obtain the following method for solving the problem (1).

Before starting calculations, we have ${\epsilon }_f>0$, $x_0\in R^n$. If $f\left(x_0\right)-f^{*}\le {\epsilon }_f$, then $x_0$ is the sought point and STOP. Otherwise, we set $h_0={\displaystyle \frac{f\left(x_0\right)-f^{*}}{\parallel \partial f\left(x_0\right)\parallel }}$, ${\xi }_0={\displaystyle \frac{\partial f\left(x_0\right)}{\parallel \partial f\left(x_0\right)\parallel }}\in R^n$, $p_0=0\in R^n$; $B_0=I$ is a unit matrix of dimension $n\times n$.

Let $x_k\in R^n$, $h_k$, ${\xi }_k\in R^n$, $p_k\in R^n$, $n\times n$ matrix $B_k$ be obtained at the iteration k.

Step 1. Calculate the next approximation

$$x_{k+1}=x_k-h_k{B}_k{\xi }_k.$$

(27)

Step 2. Calculate $f\left(x_{k+1}\right)$, $\partial f\left(x_{k+1}\right)$. If $f\left(x_{k+1}\right)-f^{*}\le {\epsilon }_f$, then $x_{k+1}$ is the sought point and STOP. Otherwise, set

$${\xi }_{k+1}={\displaystyle \frac{B_k^T\partial f\left(x_{k+1}\right)}{\parallel B_k^T\partial f\left(x_{k+1}\right)\parallel }},h_{k+1}={\displaystyle \frac{f\left(x_{k+1}\right)-f^{*}}{\parallel B_k^T\partial f\left(x_{k+1}\right)\parallel }}.$$

(28)

Step 3. Calculate

$${\lambda }_1=-{\displaystyle \frac{(p_k,{\xi }_{k+1})}{\sqrt{{(p_k,{\xi }_{k+1})}^2+{({\xi }_k,{\xi }_{k+1})}^2}}},{\lambda }_2=-{\displaystyle \frac{({\xi }_k,{\xi }_{k+1})}{\sqrt{{(p_k,{\xi }_{k+1})}^2+{({\xi }_k,{\xi }_{k+1})}^2}}}$$

and set

$$p_{k+1}=\left\{\begin{array}{ccccc}{\lambda }_1{p}_k+{\lambda }_2{\xi }_k,& \mathrm{if}& {\lambda }_1>0& \mathrm{and}& {\lambda }_2>0,\\ p_k,& \mathrm{if}& {\lambda }_1>0& \mathrm{and}& {\lambda }_2\le 0,\\ {\xi }_k,& \mathrm{if}& {\lambda }_1\le 0& \mathrm{and}& {\lambda }_2>0,\\ 0,& \mathrm{if}& {\lambda }_1\le 0& \mathrm{and}& {\lambda }_2\le 0.\end{array}\right.$$

(29)

Step 4. If $(p_{k+1},{\xi }_{k+1})\ge 0$, then $B_{k+1}=B_k$ and go to the Step 5. Otherwise, calculate

$$\eta =\left({\displaystyle \frac{1}{\sqrt{1-{(p_{k+1},{\xi }_{k+1})}^2}}}-1\right){\xi }_{k+1}-{\displaystyle \frac{(p_{k+1},{\xi }_{k+1})}{\sqrt{1-{(p_{k+1},{\xi }_{k+1})}^2}}}{p}_{k+1},$$

(30)

$$B_{k+1}=B_k\left(I+\eta {\xi }_{k+1}^T\right),h_{k+1}={\displaystyle \frac{h_{k+1}}{\sqrt{1-{(p_{k+1},{\xi }_{k+1})}^2}}},$$

(31)

$$p_{k+1}={\displaystyle \frac{1}{\sqrt{1-{(p_{k+1},{\xi }_{k+1})}^2}}}\left(p_{k+1}-(p_{k+1},{\xi }_{k+1}){\xi }_{k+1}\right).$$

(32)

Step 5. Go to the iteration $(k+1)$ with $x_{k+1}$, $B_{k+1}$, ${\xi }_{k+1}$, $h_{k+1}$, $p_{k+1}$.

The method (27)–(32) is characterized by the use of a minimum of information about the function behavior while preserving the classic Fejer step in the direction of the last calculated subgradient. In particular, in addition to the last two subgradients, it also uses an aggregate-type vector $p_k$, which is a convex combination of previously calculated subgradients, and sets a hyperplane that cuts off the set of extrema in transformed argument space from the image of the point $x_k$. In this case, the problem of updating the $p_{k+1}$ vector is solved automatically by analyzing the angles between the last subgradients and the current $p_k$.

For the method (27)–(32), as well as for the method (23)–(26), analogues of Theorems 2 and 3 hold. In particular, the following theorem is valid:

Theorem 4

([29]). Let at each step of the method (27)–(32) inequalities $\parallel B_k\parallel \le c_1$ and $\parallel \partial f\left(x_k\right)\parallel \le c_2$ hold. Then, the method (27)–(32) solves the problem (1) with an accuracy ${\epsilon }_f$ by functional with no more than in K steps, where

$$K=\left[{\left({\displaystyle \frac{c_1{c}_2\parallel x_0-x^{*}\parallel }{{\epsilon }_f}}\right)}^2\right]+1.$$

This theorem guarantees its finite convergence to the ${\epsilon }_f$ solution of the problem (1) from an arbitrary initial approximation. For ravine functions, compared to the method (23)–(26), it provides a faster tendency to zero of the matrix $B_k$ determinant as $k\to \infty$. Therefore, more effective practical work in the sense of (20) should be expected from it than from the method (23)–(26).

H form of the given methods. The methods considered above, which use one-rank ellipsoidal transformation of space, are close to methods based on two-rank ellipsoidal transformation in [27]. In particular, the first of them corresponds to the method (27)–(32) and the second to the method (37)–(44) from [27]. Let all these methods start from the same starting point. Then, thee B forms of the methods (23)–(26) and (27)–(32) and B forms of their corresponding analogues from [27] are different in the sense that they generate different sequences of matrices $B_k$. At the same time, the directions of movement in transformed spaces are also different. However, the sequences of points ${\left\{x_k\right\}}_{k=0}^{\infty }$ generated by the methods (23)–(26) and (27)–(32) coincide with the same sequences of their analogues from [27]. The fact is that the H forms of these methods coincide, that is, the sequences of ${\left\{H_k\right\}}_{k=0}^{\infty }$ matrices are the same for them. In fact, the H form of the above algorithms is determined by the expression (14) of Lemma 1 for $T_1^{-1}(\xi ,\eta ){\left(T_1^{-1}(\xi ,\eta )\right)}^T$. A similar product for $T_2(\xi ,\eta )$ has the form

$$\begin{array}{cc}\hfill & T_2^{-1}(\xi ,\eta ){\left(T_2^{-1}(\xi ,\eta )\right)}^T\hfill \\ \hfill & =R_{{\beta }_1}\left({\displaystyle \frac{\xi +\eta }{\parallel \xi +\eta \parallel }}\right)R_{{\beta }_2}\left({\displaystyle \frac{\xi -\eta }{\parallel \xi -\eta \parallel }}\right)R_{{\beta }_1}^T\left({\displaystyle \frac{\xi +\eta }{\parallel \xi +\eta \parallel }}\right)R_{{\beta }_2}^T\left({\displaystyle \frac{\xi -\eta }{\parallel \xi -\eta \parallel }}\right)\hfill \\ \hfill & =R_{{\beta }_1^2}\left({\displaystyle \frac{\xi +\eta }{\parallel \xi +\eta \parallel }}\right)R_{{\beta }_2^2}\left({\displaystyle \frac{\xi -\eta }{\parallel \xi -\eta \parallel }}\right)=\left(I+({\beta }_1^2-1)\left({\displaystyle \frac{\xi +\eta }{\parallel \xi +\eta \parallel }}\right){\left({\displaystyle \frac{\xi +\eta }{\parallel \xi +\eta \parallel }}\right)}^T\right)\hfill \\ \hfill & \times \left(I+({\beta }_2^2-1)\left({\displaystyle \frac{\xi -\eta }{\parallel \xi -\eta \parallel }}\right){\left({\displaystyle \frac{\xi -\eta }{\parallel \xi -\eta \parallel }}\right)}^T\right)\hfill \\ \hfill & =I+({\beta }_1^2-1)\left({\displaystyle \frac{\xi +\eta }{\parallel \xi +\eta \parallel }}\right){\left({\displaystyle \frac{\xi +\eta }{\parallel \xi +\eta \parallel }}\right)}^T+({\beta }_2^2-1)\left({\displaystyle \frac{\xi -\eta }{\parallel \xi -\eta \parallel }}\right){\left({\displaystyle \frac{\xi -\eta }{\parallel \xi -\eta \parallel }}\right)}^T\hfill \\ \hfill & =I-{\displaystyle \frac{(\xi ,\eta )}{2\left(1+(\xi ,\eta )\right)}}(\xi +\eta ){(\xi +\eta )}^T+{\displaystyle \frac{(\xi ,\eta )}{2\left(1-(\xi ,\eta )\right)}}(\xi -\eta ){(\xi -\eta )}^T\hfill \\ \hfill & =I+\left({\displaystyle \frac{(\xi ,\eta )}{2\left(1-(\xi ,\eta )\right)}}-{\displaystyle \frac{(\xi ,\eta )}{2\left(1+(\xi ,\eta )\right)}}\right)\left(\xi {\xi }^T+\eta {\eta }^T\right)\hfill \\ \hfill & -\left({\displaystyle \frac{(\xi ,\eta )}{2(1+(\xi ,\eta \left)\right)}}+{\displaystyle \frac{(\xi ,\eta )}{2(1-(\xi ,\eta \left)\right)}}\right)\left(\xi {\eta }^T+\eta {\xi }^T\right)\hfill \\ \hfill & =I+{\displaystyle \frac{{(\xi ,\eta )}^2}{1-{(\xi ,\eta )}^2}}\left(\xi {\xi }^T+\eta {\eta }^T\right)-{\displaystyle \frac{(\xi ,\eta )}{1-{(\xi ,\eta )}^2}}\left(\xi {\eta }^T+\eta {\xi }^T\right)\hfill \end{array}$$

and precisely corresponds to the expression (14) for $T_1^{-1}(\xi ,\eta ){\left(T_1^{-1}(\xi ,\eta )\right)}^T$.

Direct use of the relation (14) for the matrix $H_{k+1}$ recalculation requires four one-rank corrections of $n\times n$ matrices. However, by appropriately grouping the terms for one-rank matrices, two-rank corrections are enough to recalculate $H_{k+1}$. In particular, for the H form corresponding to the method (23)–(26), the formulas for the recalculation of the matrix $H_{k+1}$ take one of the forms

$$\begin{array}{cc}\hfill H_{k+1}& =H_k-{\displaystyle \frac{H_k{g}_{k+1}{g}_{k+1}^T{H}_k}{(g_{k+1},H_k{g}_{k+1})}}+{\displaystyle \frac{H_{k}p{p}^T{H}_k}{(p,H_{k}p)}},\hfill \\ \hfill p& ={\displaystyle \frac{(g_k,H_k{g}_{k+1})}{(g_k,H_k{g}_k)}}{g}_k-g_{k+1},\hfill \end{array}$$

(33)

or

$$\begin{array}{cc}\hfill H_{k+1}& =H_k-{\displaystyle \frac{H_k{g}_k{g}_k^T{H}_k}{(g_k,H_k{g}_k)}}+{\displaystyle \frac{H_{k}p{p}^T{H}_k}{(p,H_{k}p)}},\hfill \\ \hfill p& ={\displaystyle \frac{(g_k,H_k{g}_{k+1})}{(g_{k+1},H_k{g}_{k+1})}}{g}_{k+1}-g_k.\hfill \end{array}$$

(34)

Here, $g_k=\partial f\left(x_k\right)$ and $g_{k+1}=\partial f\left(x_{k+1}\right)$ are subgradients of $f\left(x\right)$ at successive points. For the H form corresponding to the method (27)–(32), the formulas for recalculating matrices $H_{k+1}$ will be a little more complicated.

In addition to saving the memory needed to store the matrix, the use of the H form also leads to more economical methods in terms of computation than the above-mentioned B methods. In particular, the two-rank correction of the matrix $H_{k+1}$, taking into account its symmetry, requires no more operations than the one-rank correction of the matrix $B_{k+1}$, and instead of two matrix–vector multiplications as required by the methods (23)–(26) and (27)–(32), one is sufficient for the H form. It is quite possible that method of justifying H methods can be no less good than B methods. However, these questions require a special study, which will not be conducted here. We only note that the formulas (33) and (34) for $H_{k+1}$ recalculation are very similar to recalculation of similar matrices in quasi-Newton methods [8].

3.3. Numerical Experiments

Since the methods (23)–(26) and (27)–(32) are identical to the methods (28)–(33) and (37)–(44) in [27], there is no particular need for a numerical study of their behavior depending on size of the problems. It can be assumed that the results given in [27] (

Table 1. ${\epsilon }_f$-convergence of the methods (23)–(26) and (27)–(32).

	Method (23)–(26)		Method (27)–(32)
Problem, n	$\mathit{iter}\left({\mathit{\epsilon }}_{\mathbf{0}}\right)$	$\mathit{iter}\left({\mathit{\epsilon }}_{\mathbf{0}}^{\mathbf{2}}\right)$	$\mathit{iter}\left({\mathit{\epsilon }}_{\mathbf{0}}\right)$	$\mathit{iter}\left({\mathit{\epsilon }}_{\mathbf{0}}^{\mathbf{2}}\right)$
$Shor$, $n=5$	112(109)	227(224)	38(36)	70(68)
$Maxquad,n=10$	120(113)	293(286)	41(35)	85(79)
$Quad\left(3\right)$, $n=5$	40(11)	73(11)	40(11)	73(11)
$Quad\left(3\right)$, $n=10$	82(60)	115(74)	76(59)	109(80)
$Quad\left(10\right)$, $n=5$	60(22)	93(22)	57(21)	90(21)
$Quad\left(10\right)$, $n=10$	187(141)	220(152)	148(124)	181(141)

and

Table 2. ${\epsilon }_f$-convergence of the method (27)–(32).

Problem, n	$\mathit{f}\left({\mathit{x}}_0\right)-{\mathit{f}}^{*}$	$\mathit{iter}\left({10}^{-5}\right)$	$\mathit{iter}\left({10}^{-10}\right)$	$\mathit{iter}\left({10}^{-20}\right)$
$Quad\left(1.1\right)$, $n=50$	581.955	42(32)	65(49)	102(73)
$Sabs\left(1.1\right)$, $n=50$	1163.909	176	279	347
$Quad\left(1.05\right)$, $n=100$	1305.013	51(43)	79(65)	124(97)
$Sabs\left(1.05\right)$, $n=100$	2610.025	318	424	614

) are valid for them. However, some specification is needed here. The fact is that the set of test problems in [27] is characterized by strong ravineness. Given that different B methods correspond to different ways of accumulating errors in the $B_k$ matrix, the numerical calculations may not completely match, that is, slight discrepancies are possible.

Therefore, in the first series of experiments, we will limit ourselves to problems of small sizes $n\sim$ 5–10 and check the numerical stability of the above methods in relation to the accuracy of solving problems by functional. As nonsmooth test problems, we will choose the two most frequently mentioned examples of piecewise-quadratic functions:

$$\begin{array}{c}Shor(n=5,f^{*}=22.6001620958)\mathrm{and}\\ Maxquad(n=10,f^{*}=-0.841408334596),\end{array}$$

which are considered poor problems for subgradient methods. Such a scrupulous indication of $f^{*}$ for the $Shor$ and $Maxquad$ problems is due to checking the methods’ work at sufficiently small values of ${\epsilon }_f$. As test smooth problems, we will choose quadratic problems $Quad\left(t\right)$ from [27] with different degrees of velocity ($Quad(10.)$ and $Quad(3.)$) and at different sizes ($n=5$ and $n=10$). Here, $f^{*}=0$. All of these problems are characterized by a single minimum point, but for nonsmooth problems, it is determined by the accuracy of the problem $f^{*}$. As ${\epsilon }_f$, we take ${\epsilon }_0$ and ${\epsilon }_0^2$. Here, ${\epsilon }_0={10}^{-5}$ for nonsmooth problems and ${10}^{-10}$ for quadratic problems.

Results of the methods work are given in Table 1. Here, $iter\left(\epsilon \right)$ is the number of calculations of the function and its subgradient values necessary to achieve $\epsilon$ accuracy by functional. The number of space transformations that are implemented in this case are given in parentheses. As can be seen from Table 1, the work of the methods (23)–(26) and (27)–(32) is not as hopeless as the work of the classic Fejer method without space transformation. This confirms to some extent that the space transformation aimed to align surface levels of the function can significantly increase the accuracy of solving problems by functional for Fejer-type methods.

The effort spent on solving these problems with the accuracy ${\epsilon }_f={\epsilon }_0^2$ by functional is often twice as much as that spent with the accuracy ${\epsilon }_f={\epsilon }_0$. However, this gap narrows significantly as the problem size increases. This is evidenced by the second series of experiments with the method (27)–(32) for weakly ravine problems $Quad\left(t\right)$ and $Sabs\left(t\right)$ from [27]. The results are shown in Table 2. For the problems $Sabs\left(t\right)$, the space transformation is implemented at almost every step of the method (27)–(32). Therefore, only the total number of iterations is given for them in Table 2.

In Figure 2 the method (27)–(32) is compared with Shor’s r-algorithm for seven unconstrained minimization problem instances with smooth and nonsmooth convex objective functions [28,41]. The number of variables varies from 10 to 50. Light grey columns in Figure 2 represent the number of iterations for r-algorithm with ${\epsilon }_x={10}^{-6}$ and ${\epsilon }_g={10}^{-6}$. Grey columns correspond to ${\epsilon }_f={10}^{-8}$ in the method (27)–(32), while dark grey columns represent results for this method with ${\epsilon }_f={10}^{-12}$. As can be seen from Figure 2, the method (27)–(32) outperforms r-algorithm in almost all seven problems except N2.

Consider two more problems to demonstrate the efficiency of Fejer-type methods. The first problem is to minimize ravine piecewise linear function $f_1(x_1,x_2)=|x_1|+t|x_2|$, $t>1$ starting from different $x_0=\left(x_0^{\left(1\right)},x_0^{\left(2\right)}\right)$. The method (27)–(32) finds the minimum point $x^{*}=(0,0)$ in no more than three iterations: 1. one iteration is required if $|x_0^{\left(2\right)}|=t|x_0^{\left(1\right)}|$, no space transformations; 2. two iterations, if $|x_0^{\left(2\right)}|<t|x_0^{\left(1\right)}|$, one space transformation; 3. three iterations if $|x_0^{\left(2\right)}|>t|x_0^{\left(1\right)}|$, one space transformation. If $|x_0^{\left(2\right)}|\ne t|x_0^{\left(1\right)}|$, then the method (2), Fejer method with Polyak’s step, converges with a geometric progression rate $q\left(t\right)=\sqrt{1-1/t^2}$ and requires a significant number of iterations if t is large. Convergence trajectories of the method (2) and the method (27)–(32) for the first problem with $t=10$ are shown in Figure 3.

The second problems is to minimize ravine piecewise quadratic function $f_2(x_1,x_2)=max\{x_1^2+{(2x_2-2)}^2-3,x_1^2+{(x_2+1)}^2\}$. The corresponding optimal solution is $x^{*}=(0,0)$ and $f_2^{*}=1$. Starting from $x_0=(1,1)$, the method (27)–(32) finds in

• 16 iterations the point $x_{16}$, where $f_2\left(x_{16}\right)-f_2^{*}\le {10}^{-6}$;
• 31 iterations the point $x_{31}$, where $f_2\left(x_{31}\right)-f_2^{*}\le {10}^{-10}$.

For the same function $f_2$, the method (2) requires more than 10,000 iterations to obtain a solution with accuracy ${\epsilon }_f={10}^{-3}$ (see Figure 4).

So, the one-rank ellipsoidal transformation of space is quite a worthy replacement of the two-rank ellipsoidal transformation in the sense that B methods based on it require fewer arithmetic operations. As with the the two-rank transformation, it allows one to build a variety of variable metric methods that have a simple geometric interpretation, and it provides a fairly convenient mechanism for proving methods in the B form based on the external approximation of the extrema set by monotonically decreasing in volume ellipsoid X, which is the image of the sphere in $Y_k$. This fact makes it possible to obtain quite satisfactory answers to a number of questions in convex programming.

In particular, Theorem 3 indicates the fundamental possibility of building a constructive stopping criterion for inconsistent systems of convex inequalities by analyzing the consequence ${\displaystyle \sum _{k=0}^{\infty }}{h}_k^2$. If ${\displaystyle \sum _{k=0}^K}{h}_k^2\ge r^2$ for some K, this means that the convex inequalities system does not have a feasible point in a sphere of radius r centered on point $x_0$. Thus, the above methods in a finite number of steps either find a feasible point of a convex inequalities system or obtain a sufficient condition for its inconsistency.

4. Orthogonalizing One-Rank Operator

The one-rank ellipsoidal operator does not allow for working with more than 2–3 vectors easily. Let $K(g_1,g_2,\dots ,g_m)$, $m\le n$ be an acute cone that localizes the extrema set. Here, $g_1$, $g_2$, …, $g_m$ are linearly independent subgradients, and $(g_i,g_j)<0$, $i,j=1,\dots ,m$, such that $i\ne j$. For effective transformation of the cone $K(g_1,g_2,\dots ,g_m)$ using the one-rank ellipsoidal operator in each of the transformed spaces, it is necessary to solve the problems of immersion of this type of cone into simpler cones $K(p_1,p_2)$, such that $(p_1,p_2)<0$. These problems are not a big trouble, and their solving requires a small computational cost. However, this can be avoided if at each step of the method we work with a simpler cone, for example, an orthogonal one. This possibility is easy to implement and the linear non-orthogonal operator described below can be used for space transformation.

Let there be a set of vectors from $R^n$

$$P=\{p_1,p_2,\dots ,p_m\},\parallel p_i\parallel =1,i=1,\dots ,m,m\le n,$$

(35)

for which the following conditions hold:

$$(p_i,p_j)=0,i\ne j,i=1,2,\dots ,m-1,j=1,2,\dots ,m-1,$$

(36)

$${∥p_m-\sum _{i=1}^{m-1}(p_m,p_i)p_i∥}^2=1-\sum _{i=1}^{m-1}{(p_m,p_i)}^2>0.$$

(37)

That is, vectors from P are linearly independent, the first $(m-1)$ vectors are mutually orthogonal, and $p_m$ is not necessarily orthogonal to all previous ones.

We denote $p={\displaystyle \sum _{i=1}^{m-1}}(p_m,p_i)p_i$. For vectors (35) satisfying (36) and (37), the following properties are fulfilled:

$$\begin{array}{ccc}\hfill \left(\mathbf{a}\right)& (p,p_m-p)=0;\hfill & \hfill \\ \hfill \left(\mathbf{b}\right)& (p_i,p_m-p)=0,i=1,\dots ,m-1;\hfill & \\ \hfill \left(\mathbf{c}\right)& (p_m,p_m-p)={\parallel p_m-p\parallel }^2.\hfill & \hfill \end{array}$$

(38)

Further, we will use them rather often.

Let us consider a linear operator from $R^n$ to $R^n$, which can be presented in matrix form as follows:

$$T_{\lambda }(p_m,p)=I+{\displaystyle \frac{p_m-p}{\parallel p_m{-p\parallel }^2}}{\left({\displaystyle \frac{1}{\lambda }}{p}_m+p\right)}^T,$$

(39)

where I is a unit matrix of dimension $n\times n$, $\lambda$ is some scalar such that $\lambda (\lambda +1)\ne 0$. We will discuss the selection of the $\lambda$ parameter and its meaning a little later. The operator (39) we will call a orthogonalizing one-rank operator. The following statements hold for it.

Lemma 4

([30]). For the operator $T_{\lambda }(p_m,p)$ for $\lambda (\lambda +1)\ne 0$, there exists an inverse $T_{\lambda }^{-1}(p_m,p)$ and

$$T_{\lambda }^{-1}(p_m,p)=I-{\displaystyle \frac{p_m-p}{\parallel p_m{-p\parallel }^2}}{\left({\displaystyle \frac{1}{\lambda +1}}{p}_m+{\displaystyle \frac{\lambda }{\lambda +1}}p\right)}^T.$$

(40)

In addition,

$$\begin{array}{cc}{T}_{\lambda }^T(p_m,p)T_{\lambda }(p_m,p)\hfill & \\ & \hfill =I+{\displaystyle \frac{1}{\parallel p_m{-p\parallel }^2}}\left(\left({\displaystyle \frac{2}{\lambda }}+{\displaystyle \frac{1}{{\lambda }^2}}\right)p_m{p}_m^T-p{p}^T+p_m{p}^T+p{p}_m^T\right).\end{array}$$

(41)

Lemma 5

([30]). Let the vector set (37) satisfy (38) and (39), and $\tilde{P}=\{{\tilde{p}}_1,{\tilde{p}}_2,\dots ,{\tilde{p}}_m\}$ be a vector set obtained from (37), as ${\tilde{p}}_i={\left(T_{\lambda }^{-1}(p_m,p)\right)}^T{p}_i$, $i=\overline{1,m}$, for $\lambda (\lambda +1)\ne 0$. Then, none of the vectors ${\tilde{p}}_i$ is identically equal to zero vector, $({\tilde{p}}_i,{\tilde{p}}_j)=0$, $i\ne j$, $i=\overline{1,m}$, $j=\overline{1,m}$, and, in addition,

$${\tilde{p}}_i=p_i,i=\overline{1,m-1},{\tilde{p}}_m={\displaystyle \frac{\lambda }{\lambda +1}}(p_m-p).$$

(42)

The proofs of Lemmas 4 and 5 are mostly based on direct verification and using the properties (38).

Therefore, the operator (39) allows for a one-rank correction of the matrix $B_{k+1}$ when moving to the next transformed space of arguments and, moreover, provides a rather convenient mechanism for work with a cone, which is defined by a finite set of vectors $m\le n$. Its use makes it possible to build constructive methods of variable metrics with a fairly visual geometric interpretation in the transformed space of arguments.

1. If one were to choose subgradients as vectors, one could obtain various methods, such as a conjugate gradients method for smooth functions. In particular, for the unconditional minimization of positive definite quadratic functions, a simple procedure of the sequential accumulation of gradients with their orthogonalization leads to methods of this type that require no more than $(n+1)$ calculations of the function and the gradient values. At the same time, using step selection, it is only necessary to guarantee the linear independence of the vectors $y_0-y_k$, $k=\overline{1,m}$, $m\le n$, where $y_0=A_m{x}_0$, $y_k=A_m{x}_k$, and, at the last step m of the method, when a linearly dependent gradient is obtained, for finding the minimum point to solve in transformed space a simple system of linear equations with a diagonal matrix of dimension $m\times m$. Controlling the $\lambda$ parameter in the operator (39) allows one to preserve positive definiteness of the $H_k=B_k{B}_k^T$ matrix and control its degeneracy degree. This allows one to build methods like the conjugate gradients method for nonsmooth problems.

2. The use of the scalar $\left(p_m,{\displaystyle \frac{p}{\parallel p\parallel }}\right)$ provides a fairly convenient mechanism for refining the lower estimate of $f^{*}$ in transformed argument space. This circumstance makes it possible to provide a constructive stopping criterion for understanding of the approximate fulfillment of the optimal conditions of convex programming problems. It can be used to construct variable metric methods using procedures like $\epsilon$-subgradient [42,43] but in a transformed argument space. Given that, at a certain value of the $\lambda$ parameter, the external approximation of the extrema set in transformed space by a sphere is preserved without increasing the radius, such procedures are also suitable for the implementation of the internal algorithm of the method of simple-body centroids [28].

3. Using the transformation (39), one can try to obtain computationally stable methods for solving systems of linear equations within the framework of the orthogonalization methods, which have proven to be numerically unstable [44,45]. At the same time, in addition to the transformation itself, there are a number of mechanisms both at the level of B form and H form for ensuring stable recalculation of the distances to hyperplanes in the transformed space of arguments.

4. The transformation (39) makes it easy to construct variable metric methods based on the localization of the extrema set by a cone in $R^n$. We will focus on one of these variants of methods for solving the problem (1), which allows for a classical Fejer step in the direction of the antisubgradient, in this section. It is based on a simple geometric fact. Let $y_k$ be the vertex of the orthogonal cone $K(g_1,\dots ,g_m)$, which localizes the extrema set, and let $g_k$ be the subgradient calculated in it. Let $K_1(g_1^{\prime },\dots ,g_{m_1}^{\prime },g_k)$ be the next cone with the vertex at $y_k$ constructed taking into account this information. Here, $g_1^{\prime },\dots ,g_{m_1}^{\prime }$ are the orthogonal subgradients of the cone $K(g_1,\dots ,g_m)$ for which the inequality $(g_i,g_k)<0$ is satisfied. By construction, this cone will contain an extrema set of the problem (1). The cone $K_1(g_1^{\prime },\dots ,g_{m_1}^{\prime },g_k)$ is non-obtuse, and its edges form acute angles. Its transformation into an orthogonal cone can naturally be used to improve the structure of the level surfaces of ravine functions.

The operator (39) and such a step-by-step cone construction strategy leads to a number of variable metric methods, which are based on the localization of the extrema set in transformed space by an orthogonal cone defined by previously calculated subgradients. At the same time, a solution of the problem of limiting the number of stored subgradients is automatically provided. Moreover, their screening is not a simple “forgetting”, but each of these vectors makes a certain contribution to space transformation, which is aimed at expanding the cone of possible directions of descent of the function.

For processes with a classical Fejer step in the direction of the antisubgradient and the described strategy of constructing a cone that localizes an extrema set, there are two possibilities for their implementation. The first is to primarily transform the space, and then, in the transformed space, to implement the classic Fejer step in the direction ${\tilde{g}}_k$. Here, ${\tilde{g}}_k$ is the image of the subgradient $g_k$ in the transformed argument space. The second is to first make a classic Fejer step in the direction $g_k$ and then transform the argument space. In both cases, we will obtain localization extrema sets by an orthogonal cone in the transformed space of arguments. In the first case, for all of the previously calculated subgradients transferred to the current point, the maximum displacement along the convexity $f\left(x\right)$ will be realized. In the second, each of the subgradients transferred to the current point, except for the last one, allows for a positive step in the direction of the antisubgradient so as not to disturb the localization of the extrema set by the orthogonal cone. So, in the second case, the method can be improved upon in the sense of the localization of the extrema set due to the step along the convex combination of orthogonal subgradients, not including the last one. This leads to a direction of movement different from the antisubgradient, and we will not consider this second option. Obviously, without a step along the convex combination of previous subgradients, it will lose to the first one.

Let us focus on the first variant and the described step-by-step strategy for building a cone. Such processes for solving the problem (1) will be called methods of orthogonal subgradient descent with a classical Fejer step due to the fact that the movement from a point in transformed space is determined by subgradient in it, which is orthogonal to the accumulated subgradients.

4.1. Orthogonal Subgradient Descent Method with a Classical Fejer Step (ORTGF)

Before proceeding to its description, we will define some additional parameters that have a purely practical meaning and specify the rule for building the next cone. In particular, in addition to the accuracy ${\epsilon }_f>0$ of solving the problem (1) by functional, we will use two more rather small positive scalars—${\epsilon }_K$ and ${\epsilon }_R$—as well as the parameter $m_0\le n-1$. Here, ${\epsilon }_K>0$ will determine the rule for screening subgradients, i.e., when moving to the $(k+1)$-th step, we will leave only those orthogonal subgradients for which the condition $\left({\displaystyle \frac{g_i}{\parallel g_i\parallel }},{\displaystyle \frac{g_k}{\parallel g_k\parallel }}\right)<-{\epsilon }_K$ is fulfilled; ${\epsilon }_R>0$ will determine additional screening of subgradients so as to avoid the accumulation of calculation errors in the orthogonalization process, that is, those orthogonal subgradients for which $\left|\left({\displaystyle \frac{g_i}{\parallel g_i\parallel }},{\displaystyle \frac{g_j}{\parallel g_j\parallel }}\right)\right|>{\epsilon }_R$ we will exclude. The $m_0$ parameter will set a limit on the maximum number of subgradients that are stored, and these will be the last computed subgradients. The set of stored subgradients at the step k will be denoted by $P_k$, and its size, that is, the number of accumulated subgradients, will be $size\left(P_k\right)$. In particular, if $P_k=\varnothing$, then $size\left(P_k\right)=0$. To work with $P_k$, we will use the sum (∪) and difference (∖) operations.

To simplify the presentation of the material, let us provide some additional restrictions in the $ORTGF$ method. In particular, we will keep only normalized subgradients and, in addition, we will consider the parameter $\lambda$ to be constant. Although, it should be noted that, by using subgradient norms and variable parameter ${\lambda }_k$, more meaningful methods of orthogonal subgradient descent can be provided. For example, by using the space dilation operator and the operator (39), it is possible to ensure the space transformation so that direction of movement in $Y_k$ coincides with the shortest vector to the convex shell of the accumulated subgradients. This direction of movement is used in $\epsilon$-subgradient methods and is considered to be a good replacement for the Newtonian direction.

Given the above, $ORTGF(x_0,{\epsilon }_f,\lambda ,{\epsilon }_K,{\epsilon }_R,m_0)$—the method of orthogonal subgradient descent with a classical Fejer step—takes the following form.

Before starting calculations, we have $B_0=I_n$, $P_0=\varnothing$.

Let $x_k$, $B_k$ and $P_k$ be obtained at the iteration k. Let us calculate $f\left(x_k\right)$ and $\partial f\left(x_k\right)$. If $f\left(x_k\right)-f^{*}\le {\epsilon }_f$, then $x_k$ is the sought point and STOP. Otherwise, we go to the iteration $(k+1)$.

Step 1. Set ${\xi }_k={\displaystyle \frac{B_k^T\partial f\left(x_k\right)}{\parallel B_k^T\partial f\left(x_k\right)\parallel }} \mathrm{and} h_k={\displaystyle \frac{f\left(x_k\right)-f^{*}}{\parallel B_k^T\partial f\left(x_k\right)\parallel }}.$

Step 2. Form a set ${\tilde{P}}_k=\{p_i\in P_k:(p_i,{\xi }_k)<-{\epsilon }_K\}$, preserving the sequence order of $p_i$ in $P_k$.

Step 3. If $size\left(P_k\right)=0$, then $B_{k+1}=B_k,{\xi }_{k+1}={\xi }_k,h_{k+1}=h_k$ and proceed to the Step 4. Otherwise, calculate the vector ${\tilde{p}}_k={\displaystyle \sum _{p_i\in {\tilde{P}}_k}}(p_i,{\xi }_k)p_i$ and recalculate the parameters $B_{k+1}=B_k{T}_{\lambda }^{-1}({\xi }_k,{\tilde{p}}_k)$,

$${\xi }_{k+1}={\displaystyle \frac{B_{k+1}^T\partial f\left(x_k\right)}{\parallel B_{k+1}^T\partial f\left(x_k\right)\parallel }},h_{k+1}={\displaystyle \frac{f\left(x_k\right)-f^{*}}{\parallel B_{k+1}^T\partial f\left(x_k\right)\parallel }}.$$

Computationally, this is equivalent to $B_{k+1}=B_k(I-{\eta }_1{\eta }_2^T)$, where

$$\begin{array}{cccc}\hfill {\eta }_1=& {\displaystyle \frac{{\xi }_k-{\tilde{p}}_k}{\parallel {\xi }_k-{\tilde{p}}_k{\parallel }^2}},\hfill & \hfill {\eta }_2=& {\displaystyle \frac{1}{\lambda +1}}{\xi }_k+{\displaystyle \frac{\lambda }{\lambda +1}}{\tilde{p}}_k.\hfill \\ \hfill {\xi }_{k+1}=& {\displaystyle \frac{{\displaystyle \frac{\lambda }{\lambda +1}}({\xi }_k-{\tilde{p}}_k)}{\parallel {\displaystyle \frac{\lambda }{\lambda +1}}(x{i}_k-{\tilde{p}}_k)\parallel }},\hfill & \hfill h_{k+1}=& {\displaystyle \frac{h_k}{\parallel {\displaystyle \frac{\lambda }{\lambda +1}}({\xi }_k-{\tilde{p}}_k)\parallel }}.\hfill \end{array}$$

Step 4. Calculate the next approximation $x_{k+1}=x_k-h_{k+1}{B}_{k+1}{\xi }_{k+1}$.

Step 5. Form the next set

$$P_{k+1}=\{p_i\in {\tilde{P}}_k:|(p_i,{\xi }_{k+1})|<{\epsilon }_R\}\cup {\xi }_{k+1},$$

preserving the sequence order of $p_i$ in ${\tilde{P}}_k$, ${\xi }_{k+1}$ will be the last vector in $P_{k+1}$. If $size\left(P_{k+1}\right)>m_0$, then $P_{k+1}=P_{k+1}\setminus \{p_1\in P_{k+1}\}$, where $p_1$ is the first vector in $P_{k+1}$.

Step 6. Go to the iteration $(k+1)$ with $x_{k+1}$, $B_{k+1}$, $P_{k+1}$.

Theorem 5

([30]). For the k-th step of the method $ORTGF$, the following inequality holds:

$$\left(A_k(x_k-x^{*}),{\xi }_k\right)\ge {\displaystyle \frac{f\left(x_k\right)-f^{*}}{\parallel B_k^T\partial f\left(x_k\right)\parallel }}\ge 0,$$

(43)

where $A_k=B_k^{-1}$. When $P_k\ne \varnothing$, then $\parallel p_i\parallel =1$ for all $p_i\in P_k$; $(p_i,p_j)=0$ for $p_i,p_j\in P_k$ such that $i\ne j$; in addition, the following inequalities hold:

$$(A_k(x_k-x^{*}),p_i)\ge 0,p_i\in P_k,$$

(44)

$$(A_k(x_k-x^{*}),{\tilde{p}}_k)\le 0.$$

(45)

To prove Theorem 5, the statements from Lemma 1 are used, as well as (44) and induction.

Theorem 6

([30]). The sequence ${\left\{x_{k+1}\right\}}_{k=0}^{\infty }$ generated by the $ORTGF$ algorithm at $\lambda =-1/2$ satisfies the inequalities

$$\parallel A_{k+1}(x_{k+1}-x^{*}){\parallel }^2\le {\parallel A_k(x_k-x^{*})\parallel }^2-{\displaystyle \frac{{(f\left(x_k\right)-f^{*})}^2}{\parallel B_k^T\partial f\left(x_k\right){\parallel }^2}}.$$

(46)

Here, $A_k=B_k^{-1}$, $A_{k+1}=B_{k+1}^{-1}$, $k=0,1,2,\dots$.

The proof of Theorem 6 is based on Theorem 5 and Lemma 4 statements.

Theorem 6 means that choosing the $\lambda =-1/2$ parameter for the $ORTGF$ method allows one to preserve the external approximation of the extrema set by a sphere with decreasing radius when moving into the transformed space of arguments. Given that $|det(B_k\left)\right|=1$, this choice of the $\lambda$ parameter allows one to substantiate the convergence of the $ORTGF$ method in the sense of the external localization of the extrema set by an ellipsoid with decreasing volume. However, the volume reduction here will be small, as it is provided only by the classic Fejer step from the center of the sphere. Nevertheless, the operator (39) at $\lambda =-1/2$ (unfortunately, it is the only one) allows for improving the structure of the function’s level surfaces, not increasing the volume of the region of localization of the extrema set. This fact makes it possible to justify variable metric methods and any other strategies for constructing a non-obtuse cone that localizes an extrema set.

In this case, the transformation (39) takes the form

$$\begin{array}{ccc}\hfill T_{-1/2}(p_m,p)& =& T_{-1/2}^{-1}(p_m,p)=I-{\displaystyle \frac{p_m-p}{\parallel p_m{-p\parallel }^2}}{\left(2p_m-p\right)}^T\hfill \\ & =& \left(I-2{\displaystyle \frac{p_m-p}{\parallel p_m-p\parallel }}{\left({\displaystyle \frac{p_m-p}{\parallel p_m-p\parallel }}\right)}^T\right)\left(I+{\displaystyle \frac{p_m-p}{\parallel p_m{-p\parallel }^2}}{p}^T\right)\hfill \end{array}$$

and contains both a “pure” mapping defined by an orthogonal matrix and a transformation acting on the “dilation” type. Relative to mapping, subgradient methods, as well as cutting-type methods, are invariant, and a transformation like “dilation” is able to improve the structure of the level surface of the function being optimized. The disadvantage of the operator (39) for $\lambda =-1/2$ is that it is not compressive for the subgradient space, i.e., $det\left(B_k{B}_k^T\right)=$1. This contributes to the accumulation of errors in the calculation of normalized directions and can lead to unstable work of methods for strongly ravine functions.

Theorem 6 implies Theorem 7.

Theorem 7

([30]). Let $\parallel B_k\parallel \le c_1$ and $\parallel \partial f\left(x_k\right)\parallel \le c_2$ at each step of the $ORTGF$ method. Then, the $ORTGF$ method solves the problem (1) with the accuracy ${\epsilon }_f$ in no more than K steps, where $K=\left[{\left({\displaystyle \frac{c_1{c}_2\parallel x_0-x^{*}\parallel }{{\epsilon }_f}}\right)}^2\right]+1$.

Its proof is identical to the proof of Theorem 3.

However, the parameter $\lambda =-1/2$ is the only one that allows, based on divergence of the sum of the series ${\displaystyle \sum _{i=0}^k}{\displaystyle \frac{{(f\left(x_i\right)-f^{*})}^2}{\parallel B_i^T\partial f\left(x_i\right){\parallel }^2}}$, one to justify convergence of the $ORTGF$ method in the sense of ${\epsilon }_f$-convergence by functional. To substantiate the convergence of $ORTGF$ for other values of the parameter $\lambda$, particularly a variable, it is necessary to use another apparatus. The interpretation of these methods as akin to conjugate directions methods [14] is most suitable here. However, this issue requires a separate discussion. Here, we limit ourselves to numerical verification of the $ORTGF$ method for two values of the $\lambda$ parameter: $\lambda =-1/2-ORTGF(-0.5)$ and $\lambda =1-ORTGF\left(1.0\right)$. However, if, for the problem (1), convergence of the $ORTGF(-0.5)$ method to the ${\epsilon }_f$ solution is justified, then the $ORTGF\left(1.0\right)$ requires a clearer justification.

The second choice of the $\lambda$ parameter can be explained by the following geometrical explanation. Let ${\tilde{P}}_k\ne \varnothing$ at the point $y_k$ at the k-th step of the $ORTGF(-0.5)$ method. Then, before transforming the cone localizing the extrema set using the operator $T_{-1/2}({\xi }_k,{\tilde{p}}_k)$, let us set ourselves the goal of reducing the volume of the region of localization of the extrema set without spoiling at the same time the structure of the cone. This is easy to do within the framework of an ellipsoid like that in [27], where one of the axes coincides with the $({\xi }_k-{\tilde{p}}_k)$ direction, and the second coincides with $p_i\in {\tilde{P}}_k$, for which $({\xi }_k,p_i)\le -1/\sqrt{2}$. This will generally be the case for strongly ravine functions. The transformation of such an ellipsoid into a sphere does not change the structure of the cone, which localizes the extrema set. It requires space dilation in the direction $({\xi }_k-{\tilde{p}}_k)$ and compression in the orthogonal direction $p_i$ and leads to a decrease of the norm of the last subgradient in a transformed argument space. Choosing the parameter $\lambda =1$ actually gives this meaning to the method, but without additional space transformations. At the same time, the norm of the last subgradient will decrease by a factor of two relative to what would be provided by the $T_{-1/2}({\xi }_k,{\tilde{p}}_k)$ operator.

4.2. Numerical Experiments

In the first series of experiments, we will check the numerical stability of the $ORTGF(-0.5)$ and $ORTGF\left(1.0\right)$ methods for the same set of test problems as in Section 2 and their ability to achieve the same accuracy ${\epsilon }_f$ by functional. As the parameters of the cone formation, we will choose: ${\epsilon }_K={10}^{-4}$, ${\epsilon }_R={10}^{-8}$, $m_0=n-1$.

The results are shown in

Table 3. ${\epsilon }_f$ convergence of the methods $ORTGF(-0.5)$ and $ORTGF\left(1.0\right)$.

	$\mathit{OTRGF}(-0.5)$		$\mathit{OTRGF}\left(1.0\right)$
Problem, n	$\mathit{iter}\left({\mathit{\epsilon }}_{\mathbf{0}}\right)$	$\mathit{iter}\left({\mathit{\epsilon }}_{\mathbf{0}}^{\mathbf{2}}\right)$	$\mathit{iter}\left({\mathit{\epsilon }}_{\mathbf{0}}\right)$	$\mathit{iter}\left({\mathit{\epsilon }}_{\mathbf{0}}^{\mathbf{2}}\right)$
$Shor$, $n=5$	33(30,4)	59(56,4)	33(30,4)	69(66,4)
$Maxquad$, $n=10$	45(37,5)	95(87,5)	42(35,5)	88(79,5)
$Quad\left(3.0\right)$, $n=5$	40(9,2)	71(9,2)	52(30,3)	96(58,3)
$Quad\left(3.0\right)$, $n=10$	80(61,5)	113(62,5)	86(68,3)	141(109,3)
$Quad\left(10.0\right)$, $n=5$	57(26,3)	90(26,3)	50(22,3)	74(22,3)
$Quad\left(10.0\right)$, $n=10$	156(123,8)	189(128,8)	131(109,4)	193(161,4)

. Here, the first number in parentheses indicates the number of transformations, and the second indicates the maximum number of accumulated subgradients. As can be seen from the table, the convergence of the methods to the ${\epsilon }_f$ solution is quite stable and almost the same as for the method (27)–(32).

Although the $ORTGF$ method is designed to align the structure of the level surface of the $f\left(x\right)$ to overcome its ravineness, it will not always achieve this goal. The point is that the ravine $\phi \left(y\right)$ at the point $y_k$ characterizes the local behavior of $\phi \left(y\right)$ near $y_k$. However, when using the classic Fejer step and the accepted way of accumulating subgradients, there is a danger that the accumulated subgradients, especially “the oldest ones”, only spoil the pattern of ravineness at the current point. For smooth functions, this will usually not work, since convergence of the method is helped by the fact that the subgradient norm tends to zero. This will ensure accumulation of subgradients that reasonably characterize the ravine. However, this is not the case for piecewise-linear functions. Here, convergence will be ensured rather by computing the vertex of the cone, which is determined by the number of pieces close to n.

Considering the above, for piecewise-linear functions, one should not expect good practical performance from $ORTGF(-0.5)$ for $m_0\ll n$. For strongly ravine functions, the situation will worsen also due to accumulation of errors in the $B_k$ matrix. Of course, the situation can be improved here both by additional space dilations aimed at volume reduction and by the procedure of “restoration” and improving the way of accumulating subgradients. For $ORTGF\left(1.0\right)$, the transformation is aimed at a stronger dilation of the cone, and it has a much better chance of good practical performance at small values of $m_0$.

In the second series of experiments, we will check the methods’ performance at minimizing a convex piecewise-linear function with a very large number of pieces. For this, we will choose “bad” for subgradient methods problem $TR48$$(n=48,f^{*}=-683,565)$ from ([46], p. 161). ${\epsilon }_K$ and ${\epsilon }_R$ are chosen the same as before. In this series of experiments, two different starting points are considered. The parameter ${\epsilon }_f=50$ for the first starting point is chosen so that the number of calculations of $f\left(x\right)$ and $\partial f\left(x\right)$ can be compared with the results presented in [47], where the points where $f\left(x\right)\le -683,500$ were considered as an acceptable solution. Results of the experiment with the first starting point show that the method $ORTGF(-0.5)$ with $m_0=n-1$ requires 139 iterations with 28 accumulated subgradients (${\epsilon }_f=50$) and 222 iterations with 31 accumulated subgradients (${\epsilon }_f={10}^{-5}$). For the second starting point, it requires 72 iterations with 34 accumulated subgradients if ${\epsilon }_f=1$ and 151 iterations with 34 accumulated subgradients if ${\epsilon }_f={10}^{-5}$. Regarding the method $ORTGF(1.)$, it was launched for values $m_0=\{20;10;5\}$, and the iteration number grows as $m_0$ decreases. In particular, the method $ORTGF(1.)$, $m_0=5$, requires 199 iterations (${\epsilon }_f=50$) and 412 iterations (${\epsilon }_f={10}^{-5}$) for the first point, and 207 iterations (${\epsilon }_f=1$) and 357 iterations (${\epsilon }_f={10}^{-5}$) for the second point, respectively. It is easy to see that the methods in the sense of ${\epsilon }_f$ convergence by functional are rather stable.

Finally, in the last, third series of experiments, we will check the numerical stability of the $ORTGF\left(1.0\right)$ method for medium-sized, strongly ravine problems ($n\sim$ 30–100). We will choose $Quad\left(t\right)$ and $Sabs\left(t\right)$ as test problems. The results are shown in

Table 4. ${\epsilon }_f$ convergence of the method $ORTGF\left(1.0\right)$ for strongly ravine problem $Quad\left(t\right)$.

Problem, n		$\mathit{Quad}\left(2.0\right)$, $\mathit{n}=30$	$\mathit{Quad}\left(1.2\right)$, $\mathit{n}=60$	$\mathit{Quad}\left(1.2\right)$, $\mathit{n}=100$
$OTRGF\left(1.0\right)$,	$iter\left({10}^{-10}\right)$	236(229,5)	188(183,5)	428(422,8)
$m_0=n-1$	$iter\left({10}^{-20}\right)$	332(315,5)	277(259,5)	542(534,8)
$OTRGF\left(1.0\right)$,	$iter\left({10}^{-10}\right)$	236(229)	188(183)	428(422)
$m_0=10$	$iter\left({10}^{-20}\right)$	332(315)	277(259)	542(534)

and

Table 5. ${\epsilon }_f$ convergence of the method $ORTGF\left(1.0\right)$ for strongly ravine problem $Sabs\left(t\right)$.

Problem, n		$\mathit{Sabs}\left(2.0\right)$, $\mathit{n}=30$	$\mathit{Sabs}\left(1.2\right)$, $\mathit{n}=60$	$\mathit{Sabs}\left(1.2\right)$, $\mathit{n}=100$
$OTRGF\left(1.0\right)$,	$iter\left({10}^{-10}\right)$	476(464,12)	464(462,21)	1480(1475,19)
$m_0=n-1$	$iter\left({10}^{-20}\right)$	527(524,12)	541(539,21)	1564(1559,19)
$OTRGF\left(1.0\right)$,	$iter\left({10}^{-10}\right)$	462(459)	469(467)	1293(1291)
$m_0=10$	$iter\left({10}^{-20}\right)$	523(520)	544(542)	1375(1373)

. Here, $iter\left(\epsilon \right)$ is the number of calculations of $f\left(x\right)$ and $\partial f\left(x\right)$. The number of transformations is given in parentheses, and for $ORTGF\left(1.0\right)$ at $m_0=n-1$, the maximum number of accumulated subgradients is given as well. If we consider n calculations of $f\left(x\right)$ and $\partial f\left(x\right)$ as one large iteration of the $ORTGF\left(1.0\right)$ method, which in terms of complexity is equivalent to one iteration of the Newton method, then, as can be seen from Table 4 and Table 5, the number of large iterations is not so big and is ∼10.

So, on the basis of the operator (39) it is also easy to construct variable metric methods with a simple geometric interpretation in a transformed argument space. At the same time, for the general problem of convex programming, it is quite possible to create practically effective methods with a strict justification of their convergence using reduction of the volume of the region of localization of the extrema set.

However, the transformation itself is only one component of the methods. The second requires choosing the direction of the function study and the method of adjusting the step multiplier. The methods used in $ORTGF$ are unlikely to be very successful, since the mandatory movement from the point of the antisubgradient greatly complicates the analysis of the $\phi \left(y\right)$ level lines. Therefore, algorithmic schemes are more rational, where the $\phi \left(y\right)$ level lines are analyzed or refined relative to a point fixed for some time. At the same time, for the calculation of $f\left(x\right)$ and $\partial f\left(x\right)$, simple procedures of the $\epsilon$-steepest descent type in the transformed space of arguments can be provided, using the orthogonality of subgradients that are preserved or their convex combinations. The issues of filtering out extra subgradients are simply solved, the upper and lower estimates of $f^{*}$ are easily specified, etc.

It is natural to direct a space transformation so that the direction of the function $f\left(x\right)$ investigation from the point approaches the direction to the optimum. Then, by changing the fixed point according to some rule (for example, to the point of record $f\left(x\right)$) or the projection of the fixed point to the set

$$Q_k=\left\{y:\phi \left(y_i\right)+\left(\partial \phi \left(y_i\right),y-y_i\right)\le f^{record},i\in I_k\right\}$$

and, repeating the same procedure with respect to a new point, we arrive at methods like Newtonian (quasi-Newton), which will be applicable to both smooth and nonsmooth functions.

5. Conclusions

Therefore, the application of the described space transformations, aimed at aligning the structure of the function level surfaces, allows us to significantly improve the performance of subgradient methods. For the convex programming problem, the use of operators (11) and (39) in combination with the space dilation operator allows us to construct a whole series of efficient variable metric methods. At the same time, for B-form methods, a simple geometric interpretation of the process in the transformed space of arguments can be given.

The considered transformations do not even come close to exhausting the list of effective one-rank space transformations. In particular, an optimal one-rank transformation like (11) is interesting, which would improve convergence in the sense of a reduction of the volume of the region of localization of an extrema set. Within the framework of orthogonalizing one-rank transformations, a simplex-like transformation is promising, which allows for building methods for convex problems of sufficiently large sizes while storing only a finite set of vectors instead of the full $n\times n$-dimensional matrix. In addition, it is possible to construct a number of other interesting linear operators of the orthogonalizing type.

Using ellipsoidal operators in Fejer-type methods results in space transformations that refine dilated level surfaces of a ravine function. This improves the behavior of these methods in further iterations even for rather large problems if the level surfaces of the objective function are not significantly dilated.

The development of nonsmooth optimization issues and applications related to solving problems of large complexity can be found in [48]. In particular, applications in economics, finance, energy saving, agriculture, biology, genetics, environmental protection, information protection, decision making, pattern recognition, and self-adaptive control of complex objects are described.

The source of high-dimensional nonsmooth problems is the problem of optimal control of systems with distributed parameters. The book [49] describes the application of nonsmooth optimization algorithms in solving some inverse problems and singular problems of optimal control of systems with distributed parameters.

Nonsmooth functions and algorithms for their minimizing are also successfully applied in such fields as optimal set partition [50,51,52], reliable networks design [53], machine learning [54], and lots of other applied areas of human practice [55,56].

Since publication of the pioneering works of N.Z. Shor, many novel applications of nondifferential optimization have been recognized in, e.g., geometrical design problems [57,58], modified Lagrangian dual problems [59,60], and green logistics [61], to mention a few. Responding to new challenges in nondifferential optimization requires novel solution techniques. Combining artificial intelligence-based approaches with classical nondifferential techniques is an interesting area for future research [62,63].