Searching for, and quantifying, non-convexity of functions
Youri Davydov Elina Moldavskaya Ričardas Zitikis









Searching for, and quantifying, non-convexity of functions





Youri Davydov^1,[youri.davydov@univ-lille1.fr],
Elina Moldavskaya^2,[elinamoldavskaya@gmail.com],
and
Ričardas Zitikis^3,[rzitikis@uwo.ca]





^1Chebyshev Laboratory, St. Petersburg State University, Vasilyevsky Island, St. Petersburg 199178, Russia



^2Department of Mathematics & Technion International School, Technion – Israel Institute of Technology, Haifa 32000, Israel




^3School of Mathematical and Statistical Sciences,
Western University, London, Ontario N6A 5B7, Canada








Abstract. Convexity plays a prominent role in a number of problems, but practical considerations frequently give rise to non-convex functions. We suggest a method for determining convex regions, and also for assessing the lack of convexity in the other regions.  The method relies on a specially constructed decomposition of symmetric matrices, such as the Hessian. We illustrate theoretical results using several examples, one of which analyses a problem arising in risk measurement and management in insurance and finance.



Key words and phrases: risk assessment, risk management, nonconvexity, penalty function, Hessian, Weyl inequality.










§ INTRODUCTION





Solutions to a great variety of optimization problems rely on convexity (e.g., Boyd and Vandenberghe, 2004), but the underlying objects (e.g., functions, surfaces, etc.) may not always be such. Several techniques can be used to overcome the challenge. One of them is based on turning non-convex objects into convex ones, and a classical example would be turning data points (e.g., incomes) into the Lorenz curve, which is convex, and has been extensively used by econometricians to measure income inequality since the pioneering works of Lorenz (1905) and Gini (1912, 1914). In turn, these ideas have given rise to lift zonoids, convex hulls, and other convex objects (e.g., Mosler, 2002). In the theory of stochastic processes, convexifications of random walks have lead to multi-dimensional convex bodies (e.g., Davydov and Vershik, 1998), as well as to convex stochastic processes (e.g., Davydov and Zitikis, 2004). While dealing with such problems, a natural question arises: how far are the obtained convexifications from the original objects?


Another technique is based on working with non-convex objects directly, perhaps initially  modifying, extending, and generalizing some of the techniques developed for tackling convex cases (e.g., Mishra, 2011). In such non-convex scenarios, the underlying functions are still convex, or concave, over certain regions of their domains of definition. A natural question arises: how can we determine, and fairly quickly due to practical considerations, those regions that are convex or concave? The present paper offers answers to questions like these by providing a rigorous methodology for determining and quantifying convexity or, in a dual way, the lack of it. The method is computationally friendly and leads to numerical assessments even when closed-form solutions are difficult to derive.


We have organized the rest of the paper as follows. In Section <ref>, we present indices of lack-of-convexity, as well as their dual versions called indices of convexity. In the same section, we provide a simple example that illustrates the concept. Section <ref> contains fundamentals that concernt the lack of positive semidefiniteness in symmetric matrices. Section <ref> is devoted to a detailed analysis of a problem that has arisen in risk measurement and management. Section <ref> concludes the paper with a brief summary of main contributions.







§ INDICES OF CONVEXITY




We start with a simple yet illuminating introduction to the main idea. Let h:(a,b)→ℝ be a real-valued function on a bounded interval (a,b) ⊂ℝ. Assume that h is  differentiable, and thus it is convex on the interval (a,b) if and only if its first derivative h' is non-decreasing on the interval. This reduces our task to the assessment of how much the derivative is increasing, and for this, we employ the idea of Davydov and Zitikis (2017). Namely, assuming that h' is differentiable, that is, the original function h is an element of the space C^2(a,b), the index of increase INC(h') of h' or, in other words, the index of convexity CONV(h) of h over the interval (a,b) is

    CONV(h)=INC(h')
    :=∫_a^b (h”(x))_+dx ∫_a^b |h”(x)| dx,

where (h”(x))_+=max{h”(x), 0}. Obviously, CONV(h)∈ [0,1] for every function h, but if h is convex on the interval (a,b), then CONV(h)=1, and if h is concave, then CONV(h)=0. This index has played a pivotal role in several applications, including financial and insurance risk management (Davydov and Zitikis, 2017), educational measurement (Chen and Zitikis, 2017), control systems assessment (Gribkova and Zitikis, 2018); we refer to Chen at al. (2018) for details and additional references on the topic. The application (details are in Section <ref>) that has inspired our present research concerns functions over multi-dimensional domains that require much more sophisticated considerations, which we describe next.

Let h:G→ℝ be a real-valued function from a bounded, open, and convex d-dimensional domain G⊂ℝ^d, for some d∈ℕ. Assume that the function is twice continuously differentiable, that is, h∈ C^2(G). Consequently, its  Hessian H_h(𝐱) exists at each point 𝐱=(x_1,… , x_d)∈ G and is a symmetric d× d matrix. The function h is convex on G if and only if the Hessian H_h(𝐱) is positive semidefinite. Hence, the problem posited in the introduction can be viewed as an assessment problem of how much, if at all, the Hessian H_h(𝐱) deviates from being positive semidefinite. In view of this, we can define the index of lack of convexity of the function h as a distance of the Hessian H_h(𝐱) from the set of all positive semidefinite, symmetric, d× d matrices. Obviously, the index is equal to 0 whenever the function h has no lack-of-convexity, that is, when the function is convex.


To proceed, we need additional notation. Let λ_1(𝐱),… , λ_d(𝐱)  denote the eigenvalues of the Hessian H_h(𝐱). They are real because the Hessian is symmetric. We define their positive and negative parts by  λ_i^+(𝐱)=max{λ_i(𝐱), 0} and λ_i^-(𝐱)=max{-λ_i(𝐱), 0}, respectively; for all i=1,… , d. Obviously, λ_i(𝐱)=λ_i^+(𝐱)-λ_i^-(𝐱) and
|λ_i(𝐱)|=λ_i^+(𝐱)+λ_i^-(𝐱).


The index of lack of convexity (LOC) of function h at point 𝐱∈ G is

    LOC(h,𝐱)=∑_i=1^d λ_i^-(𝐱) .



If the Hessian is positive semidefinite, then all its eigenvalues λ_1(𝐱),… , λ_d(𝐱)  are non-negative, and thus LOC(h,𝐱)=0, which means “zero lack of convexity.”  In other words, the function h is convex at the point x. The value of LOC(h,𝐱) never exceeds the nuclear (also known as the trace) norm ‖ H_h(𝐱) ‖_*=∑_i=1^d |λ_i(𝐱)| of the Hessian H_h(𝐱). This gives rise to our next definition.



The normalized LOC index of function h at point 𝐱∈ G is

    NLOC(h,𝐱)=∑_i=1^d λ_i^-(𝐱) ∑_i=1^d |λ_i(𝐱)|.



It follows from the definition that NLOC(h,𝐱)∈ [0,1], and thus the index is normalized. Furthermore, since NLOC(h,𝐱)=0 means no lack of convexity (i.e., convexity) and NLOC(h,𝐱)=1 means total lack of convexity, we can use the quantity
1-NLOC(h,𝐱) to measure convexity: it takes value 0 in the case of total lack of convexity and value 1 in the case of pure convexity. We have arrived at our third definition.



The index of convexity of function h at point 𝐱∈ G is

    CONV(h,𝐱)=∑_i=1^d λ_i^+(𝐱) ∑_i=1^d |λ_i(𝐱)|.





Starting from the above introduced pointwise indices of convexity, or lack of it, we can create a variety of global indices of convexity of the function h over the domain G. For example, the classical L_1-norm leads to the following global index of convexity

    CONV_1(h,G)
    =∫_G ∑_i=1^d λ_i^+(𝐱)d𝐱∫_G ∑_i=1^d |λ_i(𝐱)|d𝐱.

The index is always in the interval [0,1], takes value 0 when the function h is not convex at any point 𝐱∈ G, and takes value 1 when it is convex at every 𝐱∈ G.

Certainly, applications may suggest using other definitions of global indices, perhaps using more complex functional norms such as that of the Lebesgue space L_p(G,w,dμ) with various choices of the weight function w and positive measure μ. Nevertheless, for the sake of transparency, throughout the paper we use the L_1-norm and thus drop the subindex “1” from CONV_1(h,G) to simplify the notation and increase readability. Next is an illustrative example of the index CONV(h,G) based on a very basic yet instructive example of h.


Consider the function of two separable arguments

    h_cos(x,y)=-cos x -cos y

on the square

    S_0,0(a)=(-a,a)× (-a,a),

which is centered at (0,0) and parameterized by a∈ (0,∞ ) that we shall vary when exploring the convexity of the function h_cos over S_0,0(a). The function is depicted in Figure <ref>.

Its Hessian is the diagonal matrix whose diagonal entries are the eigenvalues λ_1(x,y)=cos x and λ_2(x,y)=cos y. Consequently, the pointwise convexity index is

    CONV(h_cos,(x,y))=(cos x )^++(cos y )^+ |cos x |+|cos y | ,

and the L_1-based global index over the square S_0,0(a) is

    CONV(a)
    :=   CONV(h_cos,S_0,0(a))
    
    
    =   ∫_-a^a ∫_-a^a ((cos x )^++(cos y )^+)dxdy ∫_-a^a ∫_-a^a (|cos x |+|cos y |)dxdy  
    
    =   ∫_0^a(cos x )^+dx ∫_0^a |cos x |dx.

The pointwise and global indices are depicted in Figure <ref>. The index CONV(a) is equal to 1 for every a ∈ (0, π/2], thus implying convexity of the function on the square S_0,0(a) for the noted a values, but when a>π/2, the amount of convexity decreases from 1 (pure convexity) to a lower, though fluctuating, level of convexity, depending on the nature of the region being absorbed by the square S_0,0(a) when the parameter a grows. This concludes Example <ref>.



In the next section we show how the aforementioned indices arise from an minimization problem in the space of matrices. The argument concerns general symmetric  matrices, and it can readily be specialized to the aforementioned Hessian.




§ INDICES OF POSITIVE SEMIDEFINITENESS





Let H be a real symmetric d× d matrix. We rewrite it as
H =Q Λ Q^⊤, where Q is an orthogonal matrix and Λ =diag(λ_1 ,… , λ_d ) is the diagonal matrix with λ_i:= λ_i(H)  denoting the eigenvalues of the matrix H. Define λ_i^+ =max{λ_i , 0} and λ_i^- =max{-λ_i , 0}. With the diagonal matrices Λ^+ :=diag(λ^+_1 ,… , λ^+_d ) and Λ^- :=diag(λ^-_1 ,… , λ^-_d ), we obtain what we call the canonical decomposition

    H =H^+ -H^-

of the matrix H, where H^+ =Q Λ^+ Q^⊤ and H^- =Q Λ^- Q^⊤. The matrices H^+ and H^- are positive semidefinite. Note also that since the nuclear norm of the matrix H is

    ‖ H  ‖_* =∑_i=1^d |λ_i |,

we have the equation ‖ H  ‖_*=‖ H^+‖_*+‖ H^-‖_*.


Of course, there can be many decompositions of H as the difference of two positive semidefinite matrices H_1 and H_2, but the canonical decomposition is a very special one, whose minimalist nature is elucidated in the following lemma.



If H=H_1-H_2 for any pair of symmetric and positive semidefinite matrices H_1 and H_2, then tr(H^+)≤tr(H_1) and also tr(H^-)≤tr(H_2).



Since H_1-H_2=H^+-H^-, the two bounds in the formulation of Lemma <ref> are equivalent. Hence, it suffices to prove only one of them, say tr(H^+)≤tr(H_1), which is equivalent to

    ∑_i=1^d λ_i^+(H) ≤∑_i=1^d λ_i(H_1)

because tr(H^+)=tr(QΛ^+Q^⊤)=tr(Λ^+).
Next we write ∑_i=1^d λ_i(H_1)=∑_i=1^d λ_i^↓(H_1), where λ_1^↓(H_1), … , λ_d^↓(H_1) denote the eigenvalues λ_1(H_1), … , λ_d(H_1) of the matrix H_1 arranged in the descending order, that is, λ_1^↓(H_1)≥…≥λ_d^↓(H_1). Next we write λ_i^↓(H_1)=λ_i^↓(H+H_2), and since H is symmetric and H_2 is positive semidefinite, we have from Weyl's (1912) monotonicity theorem (e.g., Bhatia, 1997; Corollary III.2.3, p. 63) that, for every i=1,… , d, the bound  λ_i^↓(H+H_2)≥λ_i^↓(H) holds, which is of course equivalent to the bound

    λ_i^↓(H_1)≥λ_i^↓(H).

Since H_1 is positive semidefinite, all its eigenvalues λ_1^↓(H_1), … , λ_d^↓(H_1) are non-negative. Since the matrix H is just symmetric, its eigenvalues λ_1^↓(H), … , λ_d^↓(H) may or may not be non-negative. Hence, from bound (<ref>) we have

    λ_i^↓(H_1)≥max{λ_i^↓(H),0},

which holds for all i=1,… , d. Consequently, we have

    ∑_i=1^d λ_i(H_1)=∑_i=1^d λ_i^↓(H_1)
       ≥∑_i=1^d max{λ_i^↓(H),0}
       = ∑_i=1^d max{λ_i(H),0}
    =∑_i=1^d λ_i^+(H),

which establishes bound (<ref>). The proof of Lemma <ref> is finished.



The above lemma plays a pivotal role when establishing a geometric interpretation of the index LOC(h,𝐱):  when specialized to the Hessian H(𝐱), it follows from the next theorem that LOC(h,𝐱) is the minimal nuclear-distance of the Hessian from the space ℳ^+ of all symmetric, positive semidefinite, d× d matrices.



For every  symmetric, d× d matrix H, we have

    inf_M∈ℳ^+‖ H -M ‖_* =‖ H^-‖_*.

Consequently, the following quantities

    LOPS(H)   =‖ H^-‖_*,
    
    NLOPS(H)   =‖ H^-‖_*‖ H ‖_*,
    
    PS(H)   =‖ H^+‖_*‖ H ‖_*

define, respectively, the index of lack of positive semidefiniteness (LOPS), its normalized version (NLOPS), and the index of positive semidefiniteness (PS).




We begin with the observation that since the space ℳ^+ is finite-dimensional and closed, we can find an element D∈ℳ^+ such that

    inf_M∈ℳ^+‖ H-M ‖_* = ‖ H-D ‖_*.

Denote E:=H-D and write its canonical decomposition E=E^+-E^-. Since H=D+E, we have H=D+E^+-E^- and thus H-(D+E^+)=-E^-. Due to the latter equation, we have

    ‖ H-(D+E^+) ‖_*   = ‖ E^-‖_*
       ≤‖ E ‖_*
    = inf_M∈ℳ^+‖ H-M ‖_*,

where the right-most equation holds because E=H-D. Since D+E^+∈ℳ^+, inequality (<ref>) turns into equality, and thus we have ‖ E^-‖_*= ‖ E ‖_*. This equation implies ‖ E^+‖_*=0 because ‖ E ‖_*=‖ E^+‖_*+‖ E^-‖_*, and we therefore conclude that E^+=0. Consequently, we have E=-E^-, which in turn implies -E∈ℳ^+ because E^-∈ℳ^+. Since  D∈ℳ^+ and -E∈ℳ^+, the equation H=D+E written as H=D-(-E) gives us an alternative (which may or may not be canonical) decomposition of H into the difference of two positive semidefinite, symmetric, d× d matrices.

Hence, in summary, we have two decompositions of the matrix H into the difference of positive semidefinite, symmetric, d× d matrices: the canonical decomposition H=H^+-H^- and the alternative  decomposition H=H_1-H_2 with H_1, H_2 ∈ℳ^+. We have tr(H^-)=‖ H^-‖_*. Next, upon recalling that H_2=-E, we have tr(H_2)=tr(-E)=‖ H-D ‖_*. From these observations and the bound tr(H^-)≤tr(H_2) (Lemma <ref>), we conclude that

    ‖ H^-‖_*≤‖ H-D ‖_* .

But H^- =H-H^+ and also H^+∈ℳ^+. Hence, inequality (<ref>) turns into equality, which in turn implies

    ‖ H-D ‖_* =‖ H-H^+‖_* .

Since ‖ H-H^+‖_* is equal to ‖ H^-‖_*, equation (<ref>) implies statement (<ref>) and in this way concludes the proof of Theorem <ref>.


We are now in the position to give additional insight into Definitions <ref>–<ref> by deriving alternative representations of the indices of (lack of) convexity introduced in Section <ref>.


For every h∈ C^2(G), we have the equations

    LOC(h,𝐱)=‖ H_h^-(𝐱) ‖_*,
    
    NLOC(h,𝐱)=‖ H_h^-(𝐱) ‖_*‖ H_h(𝐱) ‖_*,
    
    CONV(h,𝐱)=‖ H_h^+(𝐱) ‖_*‖ H_h(𝐱) ‖_*,

which are equivalent to equations (<ref>), (<ref>) and (<ref>), respectively.





§ AN ILLUSTRATION



The notions and quantities that we shall use in this section are standard in the literature on risk measurement and management (e.g., McNeil, Frey and Embrechts, 2015). Nevertheless, to make the text more self-contained and readable, we shall recall some of the standard definitions, such as those of the value at risk, expected shortfall, and a few other ones.




 §.§ Description



Consider a company with d business lines, i=1,… , d. Running business is costly and risky. For each i=1,… , d, let x_i be the (non-random) amount of risk capital allocated to the i^th business line in order to cover its loss X_i, which is not known beforehand and thus a random variable. Denote the cumulative distribution function (cdf) of X_i by F_i.

The expected loss not covered by the allocated capital x_i is 𝐄[(X_i-x_i)_+], where (X_i-x_i)_+ is X_i-x_i when X_i>x_i, and 0 otherwise. Obviously, if x_i is very large, then 𝐄[(X_i-x_i)_+] is very small, but keeping capital `frozen' is costly. Let ℓ_i(x_i) be the loss, or penalty, associated with the i^th business line for keeping the capital x_i. Hence, the total loss associated with the i^th business line is

    h_i(x_i)=𝐄[(X_i-x_i)_+]+ℓ_i(x_i).

It is natural to assume that ℓ_i is a non-decreasing function taking non-negative values. For example, let

    ℓ_i(x_i)=(1-p_i)x_i

for some p_i∈ (0,1). Referring to insurance, a practically relevant value would be p_i=0.99. The function h_i achieves its minimum at the point x_i^0:=VaR_p_i(X_i), where VaR_p_i(X_i) is the value-at-risk at the p_i^th quantile of the underlying loss X_i. That is,

    VaR_p_i(X_i)=inf{ x∈ℝ :  F_i(x)≥ p_i },

which, in the insurance lingo,  can be interpreted as the smallest premium that needs to be charged in order to cover at least p_i× 100% of losses. It is well known (e.g., McNeil et al., 2015) that h_i(x_i^0) is equal to (1-p_i)AVaR_p_i(X_i), where AVaR_p_i(X_i) is the average-value-at-risk at the p_i^th quantile of the underlying loss X_i, that is,

    AVaR_p_i(X_i)= 1 1-p_i∫_p_i^∞VaR_t(X_i)dt.

This is one of the fundamental risk measures currently in use in insurance and banking.


When ℓ_i(x_i) is given by equation (<ref>), the function h_i is convex, but other forms of the loss function ℓ_i are also of much interest. For example, let

    ℓ_i(x_i)=(1-p_i)x_i^α_i

for some α_i>0. Under this loss function, the function h_i may or may not be convex over its entire domain of definition, depending on the value of α_i, but it can nevertheless be convex in some regions of its domain of definition. It therefore becomes natural to specify those regions where the function h_i is convex; and where it is not, we would then wish to assess the extent of its non-convexity. The indices introduced in Section <ref> provide much needed answers and insights into such issues, and we shall illustrate this in a moment.






 §.§ Parameter choices



The losses emanating from the individual business lines need to be aggregated onto the company's level, which can successfully be implemented with the help of the weighted generalized mean

    h(𝐱)= ( ∑_i=1^d w_i h_i^β(x_i)  )^1/β

for appropriately chosen values of the parameter β∈ℝ and weights w_i≥ 0, which have to be such that ∑_i=1^d w_i=1. For example, the parameter β could be equal to -1, 0, 1 or 2, which give rise to the harmonic, geometric, arithmetic, and quadratic means, respectively. These are also the choices that we adopt in our following numerical explorations.

For the sake of expository simplicity, we deal with only two business lines, d=2, and set the two weights w_1 and w_2 to 1/2, which in practical terms means that the two business lines are viewed as being of the same importance within the company. Furthermore, let the two loss functions ℓ_i be those defined by equation (<ref>) with the parameter choices p_i=0.99, which corresponds to the 99^th risk percentile, and α_i=1/4 for both i=1,2.
Since the expected shortfall 𝐄[(X_i-x_i)_+] can be written as ∫_x_i^∞(1-F_i(t))dt, function (<ref>) becomes

    h(x,y)   = ( 1 2 ( ∫_x^∞(1-F_1(t))dt + (1-p_1)x^α_1 )^β
    +1 2 ( ∫_x^∞(1-F_2(t))dt + (1-p_2)y^α_2 )^β )^1/β
       =  ( 1 2 ( ∫_x^∞(1-F_1(t))dt + 0.01 x^1/4 )^β
    +1 2 ( ∫_x^∞(1-F_2(t))dt + 0.01 y^1/4 )^β )^1/β.

For the sake of simplicity and thus transparency, let X_1 and X_2 follow the uniform on [0,1] distribution, which means that their cdf's  F_1(t) and F_2(t) are equal to t on the unit interval [0,1]. (We can interpret X_1 and X_2 as percentages.) This turns the function h into

    h_β(x,y):= ( 1 2 g^β(x) + 1 2 g^β(y)  )^1/β

defined on the unit square [0,1]× [0,1], where
g(z)= 0.5(1-z)^2 +0.01 z^1/4 for all 0≤ z≤ 1.
We wish to explore the regions of convexity as well as of non-convexity of the function h_β, and to also quantify the extent of non-convexity when the function is not convex.


In Figures <ref> and <ref> we plot the function h_β for the aforementioned four parameter β values, that is, for β=-1, 0, 1 and 2.

Though it may look simple, the function h_β is not simple enough to easily obtain  closed-form formulas for the two eigenvalues of the Hessian, unlike we did in Example <ref> above. Hence, we resort to numerical methods, which are of primary practical interest anyway, and in addition to the surfaces themselves, in Figures <ref> and <ref> we also depict the global convexity index

    CONV(a):=CONV(h_β,S_x_0,y_0(a))
    = ∫_x_0-a^x_0+a∫_y_0-a^y_0+a(λ_1^+(x,y)+λ_2^+(x,y)) dxdy
    ∫_x_0-a^x_0+a∫_y_0-a^y_0+a( |λ_1(x,y)|+|λ_2(x,y)|) dxdy

over the square

    S_x_0,y_0(a)=(x_0,y_0)+(-a,a)× (-a,a)

and under the parameter specifications

  
  * (x_0,y_0)=(0.25,0.25) and 0<a<0.25
  
  * (x_0,y_0)=(0.25,0.75) and 0<a<0.25
  
  * (x_0,y_0)=(0.75,0.75) and 0<a<0.25
  
  * (x_0,y_0)=(0.50,0.50) and 0<a<0.5

The difference between the ranges of a in the above specifications is due to the need to keep the square
S_x_0,y_0(a) inside the domain of definition [0,1]× [0,1] of the function h_β.





 §.§ Findings



We now comment on Figures <ref> and <ref>, and start with the two  left-hand panels of Figure <ref> that correspond to the case β=-1 of the harmonic mean. Note that the three points (x_0,y_0)=(0.25,0.25), (0.5,0.5), and (0.75,0.75) are on the ridgeline of the surface, and thus result in low starting values of the global index CONV(a). Nevertheless, the index increases when the values of a increase, due to the fact that larger convex regions on both sides of the ridgeline are absorbed into the index by the expanding square S_x_0,y_0(a). Note that the functions CONV(a) corresponding to the aforementioned three points (x_0,y_0) dominate each other, which is a consequence of the fact that the points are located in places of decreasing sharpness  of the ridgeline. Note also that the function CONV(a) that corresponds to the point (x_0,y_0)=(0.25,0.75) is quite high, due to the fact that this point is far away from the ridgeline, but since the square S_x_0,y_0(a) approaches the ridgeline when a increases from 0 to 0.25, the function tends to decrease due to the absorption of increasingly larger non-convex regions.


We now look at the two right-hand panels of Figure <ref> that correspond to the case β=0.001. If compared to the previous case β=-1, the surface h_β, though still somewhat concave, it nevertheless has lost its pronounced ridgeline and thus a part of its concavity. This is reflected by higher positioned functions CONV(a) in right-hand panel (d) of Figure <ref> than those in left-hand panel (c) of the same figure. Another feature of the function CONV(a) when β=0.001 is that it drops noticeably at the end of its domain of definition, which in the case of the point (0.25,0.25) would be near a=0.25.


The functions CONV(a) depicted in panels (c) and (d) of Figure <ref> are  virtually equal to 1, except near the end-points a=0.25 and a=0.5 of the respective  domains of definition, which depend on the point (x_0,y_0). This suggests virtually total convexity of the surfaces h_β corresponding to β =1 (arithmetic mean) and β =2 (quadratic mean). Nevertheless, the functions are dropping at the end of their respective domains of definition, which suggests some loss of convexity near the edges of the unit square [0,1]× [0,1]. To understand the extent of this loss, we magnify the functions near a=0.25 and depict them in Figure <ref>.

To see where the loss of convexity happens on the surface h_β, we visualize small portions of it in Figure <ref> near the vertices (0,0) and (0,1) under the parameter choice β=1. (When β=2, the graphs are virtually the same and we therefore do not present them.) Note also from Figure <ref> that, as suggested by higher values of CONV(a) in right-hand panel (b) than in left-hand panel (a), the loss of convexity near the edges of the domain of definition [0,1]× [0,1] is smaller when β=2 than when β=1.






§ CONCLUDING NOTES



Theorem <ref> is the backbone of this paper. It introduces an index for measuring  lack of positive semidefiniteness in symmetric matrices. When specialized to the Hessian of sufficiently smooth functions, the theorem gives rise to an index suitable for measuring lack of convexity in functions, as discussed in Section <ref>. We have extensively illustrated these theoretical results in Section <ref> with the help of a non-convex function, defined on a two-dimensional convex region, that arises naturally in applications.

Since our explorations of convexity, or lack of it, rely on the Hessian and its eigenvalues, the function under consideration must have second partial derivatives. Of course, this requirement may not always be fulfilled, like for example in the two limiting cases  β↓ -∞ and β↑∞ of the function h_β defined by equation (<ref>), which give rise to non-differentiable functions min{g(x),g(y)} and max{g(x),g(y)}, respectively. To accommodate such functions, the herein developed method requires a modification, which we posit as a future problem.

It should also be noted that, from the practical point of view, the aforementioned non-differentiable limiting cases may not be necessary as one could work with functions like h_β for very small (e.g., β=-1000) or very large (e.g., β=1000) values of β, depending on the problem at hand. Though this weakening of the problem alleviates the issue associated with the existence of second partial derivatives, it comes with the necessity of employing considerable computing power.




§ ACKNOWLEDGMENTS



The research has been supported by a grant from the Natural Sciences and Engineering Research Council (NSERC) of Canada.


99

b1997
Bhatia, R. (1997).
Matrix Analysis. Springer, New York.

bv2004
Boyd, S. and Vandenberghe, L. (2004).
Convex Optimization.  Cambridge University Press, Cambridge, UK.


cz2017
Chen, L. and Zitikis, R. (2017).
Measuring and comparing student performance: a new technique for assessing  directional associations. Education Sciences, 7, 1–27.

cdgz2018
Chen, L., Davydov, Y., Gribkova, N. and Zitikis, R. (2018).
Estimating the index of increase via balancing deterministic and random data. Mathematical Methods of Statistics, 27, 83–102.


dv1998
Davydov, Yu. and Vershik, A.M. (1998).
Réarrangements convexes des marches aléatoires.
Annales de l'Institut Henri Poincare (B) Probability and Statistics, 34, 73–95.

dz2004
Davydov, Y. and Zitikis, R. (2004).
Convex rearrangements of random elements.
Fields Institute Communications, 44, 141–171.

dz2017
Davydov, Y. and Zitikis, R. (2017).
Quantifying non-monotonicity of functions and the lack of positivity in signed measures.
Modern Stochastics: Theory and Applications, 4, 219–231.


g1912
Gini, C. (1912). 
Variabilità e Mutabilità: Contributo allo Studio delle Distribuzioni e delle Relazioni Statistiche. Tipografia di Paolo Cuppini, Bologna.


g1914
Gini, C. (1914). 
On the measurement of concentration and variability of characters (English translation from Italian by Fulvio de Santis) Metron, 63, 3–38.

cdgz2018a
Gribkova, N. and Zitikis, R. (2018).
Assessing transfer functions in control systems.
Journal of Statistical Theory and Practice.
Available at 


l1905
Lorenz, M.O. (1905). Methods of measuring the concentration of wealth.
Publications of the American Statistical Association, 9, 209–219.


mfe2015
McNeil, A.J., Frey, R. and Embrechts, P. (2015).
Quantitative Risk Management: Concepts, Techniques and Tools. (Revised edition.)  Princeton University Press, Princeton, New Jersey.

m2011
Mishra, S.K. (2011).
Topics in Nonconvex Optimization: Theory and Applications.
Springer, New York.


m2002
Mosler, K. (2002).
Multivariate Dispersion, Central Regions, and Depth: The Lift Zonoid Approach.
Springer, New York.


w1992
Weyl, H. (1912).
Das asymptotische verteilungsgesetz der eigenwerte linearer partieller differentialgleichungen (mit einer anwendung auf die theorie der hohlraumstrahlung).  Mathematische Annalen, 71, 441–479.




