text:
# Learning Sparse Graph with Minimax Concave Penalty under Gaussian Markov Random Fields

This paper presents a convex-analytic framework to learn sparse graphs from data. While our problem formulation is inspired by an extension of the graphical lasso using the so-called combinatorial graph Laplacian framework, a key difference is the use of a nonconvex alternative to the ℓ_1 norm to attain graphs with better interpretability. Specifically, we use the weakly-convex minimax concave penalty (the difference between the ℓ_1 norm and the Huber function) which is known to yield sparse solutions with lower estimation bias than ℓ_1 for regression problems. In our framework, the graph Laplacian is replaced in the optimization by a linear transform of the vector corresponding to its upper triangular part. Via a reformulation relying on Moreau's decomposition, we show that overall convexity is guaranteed by introducing a quadratic function to our cost function. The problem can be solved efficiently by the primal-dual splitting method, of which the admissible conditions for provable convergence are presented. Numerical examples show that the proposed method significantly outperforms the existing graph learning methods with reasonable CPU time.

## Authors

• 1 publication
• 6 publications
• 10 publications
• 28 publications
10/22/2020

### Learning Graph Laplacian with MCP

Motivated by the observation that the ability of the ℓ_1 norm in promoti...
04/09/2019

### Sparse Laplacian Shrinkage with the Graphical Lasso Estimator for Regression Problems

This paper considers a high-dimensional linear regression problem where ...
11/30/2021

### A Convex-Nonconvex Strategy for Grouped Variable Selection

This paper deals with the grouped variable selection problem. A widely u...
10/20/2015

### Regularization vs. Relaxation: A conic optimization perspective of statistical variable selection

Variable selection is a fundamental task in statistical data analysis. S...
09/12/2018

11/06/2018

### Kernel Regression for Graph Signal Prediction in Presence of Sparse Noise

In presence of sparse noise we propose kernel regression for predicting ...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## I Introduction

How can we learn sparse graphs with enhanced interpretability under the Gaussian Markov random field (GMRF) [GMRF]? This is the central question addressed in this article. A graph, containing a set of vertices and edges, is a mathematical tool to represent the dependencies among components (such as nodes of a network or pixels of an image), through the selection of pairwise relations (edge weights) between each pair of objects (vertices). In particular, the strength of the relation can be expressed in terms of (nonnegative) graph weights. In the present context, graph “sparseness” is an important property because it tends to provide better interpretability, i.e., relative to all possible connections between nodes, only a few edges are non-zero and provide information about the major relationships between objects.

The problem of learning graphs from data has been studied widely in a variety of fields including signal processing, machine learning, and statistics

[friedman2008sparse, mazumder2012, Meinshausen_06, JMLR:v9:banerjee08a, GSP_GL_2, wang2012, Sun12, Lake, ortega16, palomar20, lam2009]. Graph learning has been considered in multiple applications such as design of functional brain network architectures [VECCHIO2017206], molecular biology [Mason2007]

, and network anomaly detection

[Bhuyan16]. We also refer the reader to [connect, Dong_servey] for comprehensive reviews of graph learning. The graphical model approach [friedman2008sparse, mazumder2012, Meinshausen_06, JMLR:v9:banerjee08a, GSP_GL_2] represents dependencies with the data in graph form and has gained significant popularity owing to two main reasons. First, the graphical model is built upon a solid statistical foundation, so that the edge weights have a physical meaning under certain assumptions. For instance, if the observed data are derived from a GMRF model, the weights are based on partial correlation coefficients [book_Probabilistic_graphical_model]. Second, it provides excellent versatility as it assumes no specific structure on the graph. A particular example of graph learning algorithm is the graphical lasso [friedman2008sparse, Bhuyan14, wang2012, Sun12, Lake], which employs regularization on the edge weights to obtain the sparse inverse covariance matrix of a GMRF model [GMRF]. This approach has been extended and modified in [ortega16] to learn several types of Laplacian matrices, including a formulation where the inverse covariance matrix has a combinatorial graph Laplacian (CGL) structure [ortega16]. As noted earlier, research on sparse graph learning is motivated by the fact that sparsity enhances the interpretability of the learned graphs [Buhlmann_book, ortega16]. All those sparsity-seeking methods exploit convex penalties (such as the norm) mainly due to their mathematical tractability.

To clarify the motivation of the present study, let us turn our attention to sparse linear regression. A plethora of nonconvex alternatives to the

regularization have been proposed to reduce the estimation bias while maintaining the benefit of variance reduction

[zhang2010, surveyOfNonconvexPenalty, Mazumder_11, lam2009]. Among them, we focus on the minimax concave (MC) penalty [zhang2010, selesnick] because:

1. it saturates (i.e., it returns a constant value when the variable being estimated exceeds a given threshold) thereby reducing the estimation biases significantly;

2. it has been shown to bridge the gap between the and norms in a parametric way [Abe_2020];

3. it is a weakly convex function [Nurminskii_73_weakly_convex]; more specifically, it is given by subtracting from the norm its Moreau envelope.

Property (iii) is of particular importance from an optimization viewpoint because the overall convexity of the cost function is ensured when the MC penalty is used with strongly convex loss functions, and also because the decomposed form in terms of a convex function and its Moreau envelope is compatible with the efficient operator splitting methods. The MC penalty has been used in various sparse estimation problems, e.g., feature selection with a sparse support vector machine (SVM)

[ex_MCP_1_Feature_Selection] and gear fault diagnosis from noisy vibration signals [ex_MCP_2_Bearing_Fault_Diagnosis]. To the best of authors’ knowledge, the quasi-norm for is the only function, excluding the MC penalty, that is known to possess property (ii) above, but it lacks properties (i) and (iii). On the other hand, the smoothly clipped absolute deviation (SCAD) penalty [lam2009] is similar to the MC penalty and it could be an alternative choice, since it actually possesses property (i) as well as the weak-convexity part of property (iii) although further investigations would be needed to determine whether property (ii) and the other part of property (iii) also hold for the SCAD penalty. While nonconvex alternatives to the penalty have been successful in the context of sparse linear regression, their study for graph learning has been limited, and most of the existing graph learning methods use the or regularization. A few exceptions include approaches using the log function [palomar20] or the norm [Shen12]. A use of the SCAD penalty was also mentioned in [lam2009, palomar_MC]. The recent works [palomar_MC, zhang2020learning] have observed that an increase of the regularization parameter of the penalty in the CGL estimation framework ultimately does not lead to a sparse solution and instead produces a dense solution associated with a fully connected graph. Based on this observation, in [palomar_MC, zhang2020learning], it has been shown that the use of the MC penalty (as well as other nonconvex penalties) yields better performance. However, these approaches are based on a nonconvex formulation and thus there is no guarantee that the generated sequence of graphs converges to a global optimum. This motivates us to devise another formulation which benefits from the weak convexity of the MC penalty to guarantee overall convexity of the entire cost so that the generated graphs converge to provably global optimum.

The goal of this article is to present a novel graph learning framework based on a convex formulation involving the nonconvex (but weakly convex) MC penalty to produce sparse graphs, and specifically sparse CGL matrices. Since CGLs are symmetric matrices, we remove this redundancy by representing a CGL matrix using a linear transform of the vector of graph weights corresponding to the upper-triangular part.111Although the one-to-one linear operator for representing the CGL was used in the structured graph learning via Laplacian spectral constraints (SGL) method [palomar20], the proposed method is more efficient (as shown in Section IV) due to the proposed reformulation, which allows to use the primal-dual splitting method [primal-dual] (as shown in Section III). Here, the upper-triangular part represents the undirected relations among nodes and completely characterizes the CGL matrix, so that our estimate is automatically guaranteed to have a Laplacian structure without the need to impose any constraints. This is in sharp contrast to the existing CGL approaches, which typically require both a positive semi-definite constraint and a linear constraint. Our formulation involves the nonconvex MC penalty, instead of the norm, while essentially keeping the same terms (the “nonsmooth” log-determinant term and the linear term) as the graphical lasso formulation but with the linear operator mentioned above. Note here that the negative log-determinant function is differentiable but with non-Lipschitz-continuous gradient. Due to the nonconvexity and the nonsmoothness, the problem cannot be solved directly using existing optimization methods. To circumvent the difficulty, we invoke the classical Moreau’s decomposition and show that the Tikhonov regularization convexifies the overall cost function, reformulating the problem into a canonical form of the primal-dual splitting method [primal-dual]. We present the admissible conditions under which the convergence to the global optimal point is guaranteed by the primal-dual splitting method. Numerical examples show that the proposed method outperforms the conventional CGL method (its -based counterpart) for three types of graph. Compared to the state-of-the-art method, the structured graph learning via Laplacian spectral constraints (SGL) [palomar20], the proposed method achieves comparable or better performance, depending on the type of graph, with up to 40 times shorter CPU time. In addition, experiments with real data show that the method produces a sparser graph than other existing methods.

New features of the present work relative to our preliminary work [koyakumaru_21] include detailed proofs of the mathematical results and refined experimental results as well as additional simulation results using real data.

## Ii Preliminaries

We present notation, and then show some mathematical tools used in this work. We finally present the primal-dual splitting method which is used to solve the proposed optimization problem to be presented in Section III.

### Ii-a Notation

The sets of real numbers and nonnegative real numbers are denoted by and , respectively. The transpose of vector/matrix is denoted by . Given a vector , define the and the norms by := and := , respectively. Similarly, given a matrix with its component denoted by , define the and the Frobenius norms by and , respectively. Given a pair of matrices and , define the inner product . Let and

denote the identity matrix and the vector of ones, respectively, and let

represent the diagonal matrix consisting of the components of a vector .

We consider undirected weighted graphs with nonnegative edge weights. The graph is composed of a set of nodes , edges , and a symmetric weight matrix with if , and if , where is the number of nodes. Here, for any by convention. CGL is defined by , where

is the degree matrix. CGL has zero row-sums with its minimum eigenvalue also zero which is simple when the graph is connected.

### Ii-B Mathematical tools

The conjugate of a function is denoted by , . The set of proper lower semicontinuous convex functions from to (] is denoted by .222 A function is proper if , and lower semicontinuous at if . The proximity operator of of index is defined as follows [combettes_book]:

 proxγf(w):=argminy∈RN(f(y)+12γ∥w−y∥22). (1)

Uniqueness and existence of the minimizer is guaranteed by the strong convexity and coercivity of . The indicator function with respect to a given set is denoted by

 ιS(w):={0,if w∈S,+∞,otherwise. (2)

It is clear by definition that . The Moreau envelope of a function of index is defined as follows [combettes_book, Definition 12.20]:

 γf(w) := miny∈RN(f(y)+12γ∥w−y∥22). (3)

Using the Moreau envelope of , which is the widely known Huber function, the MC penalty [selesnick] is defined as

 ϕMC(w) = ∥w∥1−γ∥⋅∥1(w). (4)

The nonconvex function here is known to induce a sparser and less biased estimate with respect to the penalty.

### Ii-C Primal-dual splitting method

Let and be real Hilbert spaces: in the present case, and . The primal-dual splitting method [primal-dual] solves convex optimization problems in the following form:

 minw∈X[F(w)+G(w)+H(L(w))], (5)

where is a differentiable convex function with Lipschitz continuous gradient , and are proximable proper lower semicontinuous convex functions, and is a bounded linear operator with its adjoint operator denoted by . Here, “proximable” means that the proximity operator of the function can be computed easily (in a closed form in the present case). The primal dual splitting method is given in Algorithm 1.

## Iii Proposed Algorithms

Due to its structure (i.e., symmetry and zero row-sums), the CGL is completely defined by its upper (or lower) triangular part excluding the main diagonal, or, in other words, by a length- vector, where the CGL is of size . Since all the off-diagonal components of CGL need to be nonnegative, our variable vector is constrained to the nonnegative cone (the nonnegative orthant) . Given this, we define a specific linear operator that maps a nonnegative vector of size to its corresponding CGL. For , for instance, is defined as follows:

 L:[w1,w2,w3,w4,w5,w6]T↦ ⎡⎢ ⎢ ⎢⎣w1+w2+w3−w1−w2−w3−w1w1+w4+w5−w4−w5−w2−w4w2+w4+w6−w6−w3−w5−w6w3+w5+w6⎤⎥ ⎥ ⎥⎦.

### Iii-a Problem formulation

The CGL formulation presented in [ortega16] is a popular extension of graphical lasso for imposing a Laplacian constraint. With a slight modification using the linear operator introduced above, the CGL formulation is given by

 P0:minw∈C−logdet(L(w)+J)+⟨S,L(w)⟩+λ1∥w∥1, (6)

where , stands for the sample covariance obtained from data, and is the regularization parameter. Note here that is positive definite if and only if the graph of is connected, i.e., is an irreducible matrix, where denotes the Hadamard product.

Replacing the regularization term of Problem by the MC penalty given in (4), the problem reads as follows:333Although the formulation in P has been considered in the literature [palomar_MC, zhang2020learning], it was also considered earlier in the authors’ previous works [koyakumaru_bthesis, koyakumaru_21] as an intermediate step, and the present study is independent from [palomar_MC, zhang2020learning].

 P1:minw∈C−logdet(L(w)+J)+⟨S,L(w)⟩ (7) +λ1[∥w∥1−γ∥⋅∥1(w)].MC

We introduce the Tikhonov regularization term , , which also plays a role of convexification as shown below. Introducing the indicator function to accommodate the constraint as well, Problem is transformed into the following unconstrained optimization problem:

 P2: minw∈Rn(n−1)2−λ1 γ∥⋅∥1(w)+λ22∥w∥22F(w) +ιC(w)+λ1∥w∥1+⟨S,L(w)⟩G(w)−logdet(L(w)+J)H(L(w)). (8)

By using Moreau’s decomposition [combettes_book, Theorem 14.3], the MC penalty term can be rewritten as

 −γ∥⋅∥1(w) = γ−1(∥⋅∥∗1)(γ−1w)−12γ∥w∥22 (9) = γιB∞(w)−12γ∥w∥22, (10)

where is the ball of radius- [combettes_book, Example 13.32]. The usefulness of this decomposition of the MC penalty term has been observed also in [suzuki21, kaneko20, yukawa21_alime, komuro_2021_SSP]. Using (10), the function of P can be rewritten as

 F(w)=λ1 γιB∞(w)−λ12γ∥w∥22+λ22∥w∥22, (11)

of which the convexity is ensured clearly by choosing and such that (see Proposition 4 below). The Tikhonov regularization term thus has a convexification property, as mentioned above. On the other hand, the functions and are convex, since the composition of a convex function with an arbitrary affine operator is also a convex function. Hence, under the convexity condition given above, Problem takes the form of (5), and it can be solved by the primal-dual splitting method.

### Iii-B Optimization algorithm

The proposed algorithm is derived by applying the primal-dual spitting method to Problem .

#### Iii-B1 Derivation of ~wk+1

We define the soft thresholding operator for a length- positive vector by

 [softδ(w)]i=⎧⎪⎨⎪⎩wi−δi,if wi≥δi,0,if |wi|<δi,wi+δi,if wi≤−δi, (12)

where is the th component of the argument. The convex projection onto the nonnegative cone is given by

 [PC(w)]i={wi,if wi≥0,0,if wi<0. (13)

Applying Step 1 of Algorithm 1 to Problem yields

 ~wk+1 =proxτG[wk−τL∗(Vk)−τ∇F(wk)] (14) =proxτG[wk−τL∗(Vk)−τ(γ−1λ1prox∥⋅∥1(wk) (15)

The operators and can be computed by using the following propositions.

###### Proposition 1

Let be an arbitrary CGL matrix with its component denoted by . Then, for any such that , it holds that

 [L∗(M)](2n−p−1)p/2+q−n=mp,p+mq,q−mp,q−mq,p.

Proof: See Appendix A.

###### Proposition 2

The proximity operator of of index can be expressed by

 proxτG(w)=PC(w−τ(λ11+L∗(S))). (16)

Proof: See Appendix B.

By using Proposition 2 and , (15) can be rewritten as

 ~wk+1=PC[wk−τL∗(Vk)−τ(λ11+L∗(S)) −τ(γ−1λ1soft1(wk)−γ−1λ1wk+λ2wk)]. (17)

#### Iii-B2 Derivation of ~Vk+1

Substituting Moreau’s decomposition [combettes_book, Theorem 14.3]

 proxσH∗(u) = u−σproxσ−1H(σ−1u) (18)

with into Step 2 of Algorithm 1 yields

 ~Vk+1 =Vk+σL(2~wk+1−wk) − σproxσ−1(−logdet(⋅+J))[σ−1Vk+L(2~wk+1−wk)].

Here, the proximity operator can be written in a closed form, as shown in the following proposition.

###### Proposition 3

For a positive semi-definite matrix , it holds that

 proxσ−1(−logdet(⋅+J))(W) =Qdiag⎛⎜ ⎜⎝μ1+√μ21+4σ−12,⋯,μn+√μ2n+4σ−12⎞⎟ ⎟⎠QT −J, (20)

where is the th eigenvalue of , and

with the eigenvectors

of .
Proof: See Appendix C.

An application of Proposition 3 to (LABEL:y_np1_op) yields

 ~Vk+1=Vk+σL(2~wk+1−wk)+σJ −σ⎡⎢ ⎢⎣Udiag⎛⎜ ⎜⎝ν1+√ν21+4σ−12,⋯,νn+√ν2n+4σ−12⎞⎟ ⎟⎠UT⎤⎥ ⎥⎦, (21)

where is the eigenvalue of , and with the eigenvectors of .
The proposed algorithm is given in Algorithm 2. Our formulation based on the MC penalty is expected to yield a sparser solution with better interpretability than the conventional -based methods due to the efficient sparsity promoting property of the MC penalty. In addition, our representation of CGL using the linear operator reduces the number of variables approximately by half, while also transforming the positive semi-definite constraint of the graph Laplacian to the nonnegativity constraint . However, our formulation needs complexity due to the need to execute matrix multiplication and eigenvalue decomposition. This computational drawback can be mitigated by using the eigenvalue decomposition method for symmetric matrices [MagoarouGT16]. Note that, in the particular case of , the proposed algorithm gives an alternative way to solve the graphical lasso problem for CGL.

### Iii-C Convergence conditions

Convergence is guaranteed under the following conditions.

###### Proposition 4

If , the function is convex. In this case, Algorithm 1 converges to a global minimizer of if the following conditions are jointly satisfied:

Proof: See Appendix D.

From a theoretical side, a use of satisfying the convexity condition shown in Proposition 4 ensures convergence to a global minimizer. From a practical side, on the other hand, a use of violating the convexity condition may yield better performance, as will be seen in Section IV. However, we emphasize that improved performance comes with no theoretical guarantees. A remarkable advantage of the present framework is its flexibility due to the use of powerful convex analytic solver, which allows to extend the presented framework in many possible directions including dynamic graph learning [HallacPBL17, tanaka19_time_varying].

## Iv Numerical Examples

We show the efficacy of the proposed method through some experiments with synthetic and real data. We first show the performances of the proposed method for different regularization parameters. We then compare the performance of the proposed method with CGL estimation [ortega16] and SGL estimation [palomar20]. 444We used the implementations of the SGL and CGL algorithms in spectralGraphTopology (https://CRAN.R-project.org/package=spectralGraphTopology).

### Iv-a Experiments with synthetic data

Dataset generation: We consider three types of graph: (i) grid graph with nodes connected to their four nearest neighbors (except the nodes at boundaries), (ii) random modular graph (a.k.a. stochastic block model)

with four modules where the nodes are connected across the modules and within each module with probabilities 0.01 and 0.3, respectively, and (iii) Erdös-Rényi graph

with nodes connected to other nodes with probability 0.1. The graph weights are randomly drawn from the uniform distribution over the interval [0.1, 3.0], regarded as the ground-truth graph Laplacian

in this experiment. From each graph generated, data are generated from , where denotes the Moore-Penrose pseudo inverse, and the covariance matrix is computed from data. For each type of graph, we randomly generate graphs with nodes using the toolbox given in [perraudin2014gspbox].

Performance measure:

The relative error (RE) and F-score (FS) are used as performance measures:

 RE(ˆΘ,Θ⋆):=∥∥ˆΘ−Θ⋆∥∥2F∥Θ⋆∥2F, (22) FS(ˆΘ,Θ⋆):=2tp2tp+fn+fp, (23)

where , and stand for true-positive, false-positive, and false-negative, respectively. Here, RE indicates the discrepancy between the ground-truth graph Laplacian and its estimate , while FS is a measure of accuracy of binary classification (taking values in [0,1]), indicating whether the sparse structures are extracted correctly.

#### Iv-A1 Performance of the proposed method

We study the impacts of the parameters and of the MC penalty on the performance of the proposed method. We also tested the case of , which makes the entire cost function nonconvex for any , with the other parameters tuned manually (see Table I). To study the impact of , we fix which gave a best performance in the nonconvex case. Since is the algorithm parameter and it only affects the convergence speed, we fix it to . The other parameters are then set to and according to the convexity condition (See Proposition 4). Figure 1 plots the RE and FS curves across in modular graph for different choices of . To study the impact of , on the other hand, we fix and choose the other parameters in the same way as in Fig. 1. Figure 2 plots the RE and FS curves for different choices of under the same conditions as in Fig. 1, with =. In Figs. 1 and 2, the proposed method using the convexity condition attains better performance than the nonconvex case when is small, while the nonconvex case is better when is large. To be more specific, when is small, using larger , or larger , yields better performance. Although the regularization parameter for the Tikhonov regularization needs to be sufficiently large to ensure the convexity of the entire objective, using a that is too large tends to yield a less sparse solution, which means degradation of graph interpretability (cf. Section III-C).

#### Iv-A2 Comparisons in estimation accuracy

Parameters: The best parameters are chosen manually, see Table I. The generated graphs depend only on , , and , while being independent of the algorithm parameters , , and in principle. It is very important to tune the parameters , , and carefully for better performance, although the generated graph changes gradually as each of those parameters changes. As shown in Section IV-A1, for the proposed method under the convexity condition, is used. We mention that a smaller value of tends to give a better result for RE and FS when is large. Although the convergence is guaranteed under the condition of described in Section III-C, (for which there is no guarantee of convergence to the global minimizer) gave the best performance in the current experiments. Regarding the parameters for CGL and SGL, we follow the parameter selection techniques proposed in [ortega16] and [palomar20], respectively.

Results: Figure 3 shows the ground truth and the learned graph for where is the number of measurements. One can see that the proposed method yields a more accurate graph than the other methods; in particular, the graph obtained by the proposed method is remarkably sparse. Figures 46 show the performances in RE and FS across for the grid graph , the random modular graph , and the Erdös-Rényi graph , respectively. In Fig. 4, the proposed method (nonconvex) significantly outperforms CGL due to the use of the MC penalty, while the proposed method under the convexity condition achieves approximately the same RE performance as the nonconvex case with the degraded FS performance for large . The performance of SGL is close to that of the proposed method in this case. The proposed method under the convexity condition exhibits approximately the same performance as the proposed method for (the nonconvex case). In Fig. 5, on the other hand, the proposed algorithm outperforms SGL considerably. The proposed method significantly outperforms CGL in FS for large , although those two methods exhibit comparable performances in RE (and in FS as well for small ). In Fig. 6, the proposed method outperforms the other methods in FS for a wide range of values, while the proposed method under the convexity condition achieves approximately the same RE performance as the nonconvex case with its FS performance close to those of CGL and SGL for large . We remark that the difference between the proposed method and CGL is more notable in Fig. 4 than in Figs. 5 and 6 because the graph in Fig. 4 is approximately four times sparser than that of Fig. 5 and twice sparser than that of Fig. 6. Thus, the sparsity assumption is a better match to the actual data in for the graph from Fig. 4. We finally remark that, for CGL and SGL, a small regularization parameter was used because use of large regularization parameters with the norm leads to increased errors, which degrade the quality (e.g., interpretability) of the learned graphs. This is the reason why the graphs obtained by CGL and SGL in Fig. 3 are not sufficiently sparse.

#### Iv-A3 Comparisons in computation time

We investigate how the computation time with tolerance error changes with the size of the graph for the modular graph We set and perform graph learning 15 times as in Section IV-A. The computation time for is summarized in Table II. It can be seen that the proposed method is 5.38 – 43.4 times faster than SGL. Although the computation time of the proposed method is higher due mainly to the eigenvalue decomposition of the larger sized matrix compared to CGL, we emphasize that the performance improvements are remarkable especially for the grid graph (See Fig. 4). The significant advantage in CPU time is due to the fact that the proposed method requires a few thousand iterations for approximate convergence, while SGL requires over iterations on average with per-iteration complexity of order .

### Iv-B Experiments with real data

We test our method for the animal dataset [Kemp10687], in which each node represents each animal and the edges represent how much the animals are related to each other. The dataset contains binary values (i.e., it is a categorical dataset) which represent the answers to some questions such as “has lungs?” for instance. There are 102 such questions in total, answered for 33 different animals. The covariance matrix is created based on this data set, and the graph is learned in the same way as in the previous experiment. The results are shown in Fig. 7. It can be seen that the proposed method produces a sparser graph than CGL and SGL while preserving the dominant links.

## V Conclusion

We presented a graph learning method inserting the nonconvex MC penalty into the extension of the graphical lasso formulation to produce sparse and accurate graphs. With the linear-operator-based representation of CGL together with the reformulation through the Moreau decomposition, an efficient algorithm was derived with the primal-dual splitting method, for which an admissible choice of parameters was presented to ensure the provable convergence. Numerical examples showed that the proposed method significantly outperformed CGL for the high-sparsity grid graphs, while still achieving higher F-scores than CGL and SGL (the state-of-the-art method) for the low-sparsity modular graphs. In addition, the comparisons in CPU time showed that the proposed method was dramatically faster than SGL.

## Appendix A Proof of Proposition 1

Let . By definition of adjoint operator, it holds that =, of which the left side can be expanded as . It can therefore be seen that and . In general, it can be verified that .

## Appendix B Proof of Proposition 2

The assertion can be verified by combining the basic property [combettes_book, Proposition 24.8]

 proxϕ+⟨⋅,u⟩(x)=proxϕ(x−u) (B.1 )

and the following fact [yukawa, Proposition 1]:

 proxιC+η∥⋅∥1(w)=PC(w−η1) (B.2 )

for . Indeed, since , we have . Applying (B.1) to and , we obtain

 proxτG(w) = proxιC+τλ1∥⋅∥1+⟨⋅,τL∗(S)⟩(w) (B.3 ) = proxιC+τλ1∥⋅∥1(w−τL∗(S)) = PC(w−τL∗(S)−τλ11).

## Appendix C Proof of Proposition 3

The proximity operator of

 f(W)↦{−logdet(W), if W≻0,+∞, otherwise,

is given by [combettes_book, Example 24.66]

 prox−σ−1logdet(⋅)(W)=12Qdiag(μ1+√μ21+4σ−1,⋯,μn+√μ2n+4σ−1)QT.

By (C) and the property [combettes_book, Proposition 24.8]:

 proxf(⋅+J)(W)=proxf(W+J)−J, (C.2 )

we obtain the result.

## Appendix D Proof of Proposition 4

It is clear from (11) that is convex when . The convergence condition of the primal dual splitting method is given as follows [primal-dual]:

1. ,

2. ,

where is the Lipschitz constant of and is the operator norm. We show below the Lipschitz constant of and the operator norm of . (Although it is shown in [palomar20] that , we show the proof for self-containedness.)

### D-a Derivation of the Lipschitz constant of ∇ F

The gradient of is given by

 ∇F(w)=γ−1λ1prox∥⋅∥1(w)+(λ2−γ−1λ1)w.

Hence, by the nonexpansiity of the proximity operator as well as the triangular inequality of norm, we obtain

 ∥∇F(w)−∇F(w′)∥2 =∥∥(γ−1λ1prox∥⋅∥1(w)+(λ2−γ−1λ1)w) −(γ−1λ1prox∥⋅∥1(w′)+(λ2−γ−1λ1)w′)∥∥2 (D.1 ) ≤∥γ−1λ1prox∥⋅∥1(w)−γ−1λ1prox∥⋅∥1(w′)∥2 +∥(λ2−γ−1λ1)(w−w′)∥2 (D.2 ) ≤γ−1λ1∥w−w′∥2+(λ2−γ−1λ1)∥w−w′∥2 (D.3 ) ≤λ2∥w−w′∥2, (D.4 )

from which is -Lipschitz contiuous.

### D-B Derivation of ∥L∥

By definition of Laplacian, we have


text:
# Tips and tricks¶

## Making models fit on Loihi¶

### Splitting large Ensembles¶

By default, NengoLoihi will split Ensemble objects that are too large to fit on a single Loihi core into smaller pieces to distribute across multiple cores. For some networks (e.g. most densely-connected networks), this can happen by itself without any guidance from the user.

For networks that use nengo.Convolution transforms, such as image processing networks, some assistance is usually required to tell NengoLoihi how to split an ensemble. This is because grouping the neurons sequentially is rarely ideal for such networks. For example, if an Ensemble is representing a 32 x 32 x 4 image (that is 32 rows, 32 columns, and 4 channels), we might want to split that ensemble into four 32 x 32 x 1 groups, or four 16 x 16 x 4 groups. In the first case, each Loihi core will contain information from all spatial locations in the image, but each will only contain one of the channels. In the second case, each Loihi core will represent a different spacial quadrant of the image (i.e., top-left, top-right, bottom-left, bottom-right), but each will contain all channels for its respective location. In neither case, though, will all cores contain solely consecutive pixels, assuming our pixel array is ordered by rows, then columns, then channels. Since the default behaviour of the system is to split into consecutive groups, we need to override that behaviour.

To do this, we use the BlockShape class with the block_shape configuration option:

import numpy as np

image_shape = (32, 32, 4)

with nengo.Network() as net:

# first case: splitting across channels
ens1 = nengo.Ensemble(np.prod(image_shape), 1)
net.config[ens1].block_shape = nengo_loihi.BlockShape((32, 32, 1), image_shape)

# second case: splitting spatially
ens2 = nengo.Ensemble(np.prod(image_shape), 1)
net.config[ens2].block_shape = nengo_loihi.BlockShape((16, 16, 4), image_shape)

We are not limited to splitting only along the spatial or channel axes. For example, with a 32 x 32 x 4 image we could choose a block shape of 16 x 16 x 2, which would result in 8 cores tiling the image both in the spatial and channel dimensions. We could also use shapes that are uneven in the spatial dimensions, for example 16 x 32 x 2.

Furthermore, the block shape does not have to fit evenly into the image shape in all (or even any) of the dimensions. For example, with a 4 x 4 image, we could choose a block shape of 3 x 3; this would result in 4 blocks: a 3 x 3 block for the top-left of the image, a 3 x 1 block for the top-right, a 1 x 3 block for the bottom-left, and a 1 x 1 block for the bottom-right. In this case, it would be better to use a 2 x 2 block shape, which also results in 4 blocks, but uses resources more equally across all cores. (This assumes that resource constraints are preventing us from using e.g. a 2 x 4 or 4 x 4 block shape that would simply use fewer cores.)

The constraints on BlockShape are that each block has to fit on one Loihi core. The most basic resource limitation is that the number of neurons (the product of the shape), must be less than or equal to 1024 (the maximum number of neurons per core). Our two original block shapes of 32 x 32 x 1 and 16 x 16 x 4 both equal exactly 1024 neurons per core. However, there are other limiting resources on Loihi cores, such as the numbers of input and output axons, and the amount of synapse memory. We therefore may not always be able to use block shapes that fully utilize the number of compartments, if other resources are in short supply.

### Measuring utilization of chip resources¶

The Model.utilization_summary command can be used to get more information on the resources used by each block (i.e. Loihi core). This can help you to judge whether cores are being optimally utilized. When using this command, it is best to give all your ensembles unique labels (like nengo.Ensemble(..., label="my_ensemble")); these names will show up in the summary, allowing you to identify problematic blocks.

### Reducing axons by changing pop_type¶

The pop_type configuration option can be used to set the type of population axons used on convolutional connections. Setting this to 16 instead of 32 (the default) reduces the number of axons required by the model, but also adds some restrictions on how convolutional connections are set up. See add_params for more details.

### Reducing synapses and axons by changing block shape¶

In networks with convolutional connections, an inefficient parameterization can cause some connection weights to be copied up to four times (to work around limitations when mapping these connections onto Loihi). If you are running out of synapse memory or axons on convolutional connections, use the following guidelines to see whether restructuring could help.

First, it is important to understand the difference between the spatial dimensions of a shape and the channel dimension. The channel dimension will be the first dimension of the shape if channels_last=False or the last dimension if channels_last=True. All other dimensions of the shape will be part of the spatial shape. For example, if channels_last=True and we have the shape 32 x 32 x 4, the spatial shape is 32 x 32, and the spatial size is 32 * 32 = 1024.

When choosing our block shape, we can trade off between spatial size and channel size. For example, if our image shape is 32 x 32 x 8, an Ensemble representing this will need at least 8 cores (since 32 * 32 * 8 / 1024 == 8). Two potential block shapes to achieve these 8 cores are 32 x 32 x 1, which has a spatial size of 1024, and 16 x 16 x 4, which has a spatial size of 256. If we wish to reduce the spatial size of our block shape, decreasing the spatial size per block while simultaneously increasing the number of channels per block will often let us keep the same number of neurons per core, but decrease other resources such as synapse memory or axon usage.

If the problem is output axon usage, try to increase the channel size of the blocks on any Ensemble targeted by connections from the problematic Ensemble. Axons can be reused across channels, so the more channels per core, the fewer output axons are required to send the information to all target cores.

If the problem is input axon usage, try to reduce the spatial size of the problematic Ensemble. Again, since axons are reused across channels, changing the channel size will have no effect (potentially, it can even be increased to compensate for the drop in spatial size, and keep the same number of compartments per core).

If the problem is synapse memory usage, then the problem is caused by incoming Connections to the problematic Ensemble. The solution depends on the value of channels_last on the Convolution transform, and the value of pop_type on the Connection (if you have not set pop_type, 32 is the default value). The following can be applied to any or all of the incoming Connections:

• If the connection is using channels_last=False and pop_type=32, extra weights are created if the spatial size is greater than 256 (the factor by which the size of the weights is multiplied is approximately the spatial size divided by 256). Decrease the spatial size.

• If the connection is using channels_last=False and pop_type=16, extra weights are always created. Consider using channels_last=True, or not using pop_type=16 if you are using less than 50% of the available axons.

• If the connection is using channels_last=True and pop_type=32, extra weights are created if there are more than 256 neurons per core. Consider using channels_last=False.

• If the connection is using channels_last=True and pop_type=16, extra weights are created if the number of channels per block is not a multiple of 4, and if there are more than 256 neurons per core. Consider making the channels per block a multiple of 4.

In all cases, decreasing the number of channels per block will decrease the amount of synapse memory used, since there is one set of weights per channel.

## Local machine¶

### SSH hosts¶

Adding ssh hosts to your SSH configuration will make working with remote superhosts, hosts, and boards much quicker and easier. After setting them up, you will be able to connect to any machine through a single ssh <machine> command.

To begin, make a ~/.ssh/config file.

touch ~/.ssh/config

Then open that file in a text editor and add a Host entry for each machine that you want to interact with remotely.

Typically machines that you can connect to directly will have a configuration like this:

Host <short name>
HostName <host name or IP address>

For security, the port on which ssh connections are accepted is often changed. To specify a port, add the following to the Host entry.

Host <short name>
...
Port 1234

Finally, many machines (especially hosts and boards) are not accessible through the open internet and must instead be accessed through another machine, like a superhost. To access these with one command, add the following to the Host entry. <tunnel short name> refers to the <short name> of the Host entry through which you access the machine (e.g., the <host short name> entry uses the superhost’s short name for <tunnel short name>).

Host <short name>
...
ProxyCommand ssh <tunnel short name> -W %h:%p

Once host entries are defined, you can access those machine with:

ssh <short name>

You can also use the short name in rsync, scp, and other commands that use ssh under the hood.

For more details and options, see this tutorial.

We recommend that Loihi system administrators make specific host entries for their system available to all users.

### SSH keys¶

SSH keys allow you to log in to remote machines without providing your password. This is especially useful when accessing a board through a host and superhost, each of which require authentication.

You may already have created an SSH key for another purpose. By default, SSH keys are stored as

• ~/.ssh/id_rsa (private key)

• ~/.ssh/id_rsa.pub (public key)

If these files exist when you do ls ~/.ssh, then you already have an SSH key.

If you do not have an SSH key, you can create one with

ssh-keygen

Follow the prompts, using the default values when unsure. We recommend setting a passphrase in case someone obtains your SSH key pair.

Once you have an SSH key pair, you will copy your public key to each machine you want to log into without a password.

ssh-copy-id <host short name>

<host short name> is the name you specified in your SSH config file for that host (e.g., ssh-copy-id loihi-host). You will be prompted for your password in order to copy the key. Once it is copied, try ssh <host short name> to confirm that you can log in without providing a password.

### Remote port tunneling¶

Tunneling a remote port to your local machine allows you to run the Jupyter notebook server or the NengoGUI server on the superhost or host, but access the web-based interface on your local machine.

To do this, we will create a new terminal window on the local machine that we will keep open while the tunnel is active. In this terminal, do

ssh -L <local port>:localhost:<remote port>

You will then enter an SSH session in which you can start the process that will communicate over <remote port>.

Example 1: Starting a NengoGUI server on port 8000 of superhost-1, which has a loihi conda environment.

# In a new terminal window on your local machine
ssh -L 8000:localhost:8000 superhost-1
# We are now on superhost-1
source activate loihi
cd ~/nengo-loihi/docs/examples
nengo --port 8000 --no-browser --auto-shutdown 0 --backend nengo_loihi

On your local machine, open http://localhost:8000/ and you should see the NengoGUI interface.

Example 2: Starting a Jupyter notebook server on port 8080 of superhost-2, which has a loihi virtualenv environment.

# In a new terminal window on your local machine
ssh -L 8080:localhost:8080 superhost-2
# We are now on superhost-2
workon loihi
cd ~/nengo-loihi/docs/examples
jupyter notebook --no-browser --port 8080

The jupyter command should print out a URL of the form http://localhost:8888/?token=<long-strong>, which you can open on your local machine.

### Syncing with rsync¶

If you work on your local machine and push changes to multiple remote superhosts, it is worth spending some time to set up a robust solution for syncing files between your local machine and the superhosts.

rsync is a good option because it is fast (it detects what has changed and only sends changes) and can be configured to ensure that the files on your local machine are the canonical files and are not overwritten by changes made on remotes. rsync also uses SSH under the hood, so the SSH hosts you set up previously can be used.

rsync is available from most package managers (e.g. apt, brew) and in many cases will already be installed on your system.

The basic command that is most useful is

rsync -rtuv --exclude=*.pyc /src/folder /dst/folder

• -r recurses into subdirectories

• -t copies and updates file modifications times

• -u replaces files with the most up-to-date version as determined by modification time

• -v adds more console output to see what has changed

• --exclude=*.pyc ensures that *.pyc files are not copied

When sending files to a remote host, you may also want to use the --delete option to delete files in the destination folder that have been removed from the source folder.

To simplify rsync usage, you can make small bash functions to make your workflow explicit.

For example, the following bash functions will sync the NxSDK and nengo-loihi folders between the local machine and the user’s home directory on host-1. In this example, the --delete flag is only used on pushing so that files are never deleted from the local machine. The --exclude=*.pyc flag is only used for nengo-loihi because *.pyc files are an important part of the NxSDK source tree. These and other options can be adapted based on your personal workflow.

LOIHI="/path/to/nengo-loihi/"
NXSDK="/path/to/NxSDK/"
push_host1() {
rsync -rtuv --exclude=*.pyc --delete "$LOIHI" "host-1:nengo-loihi" rsync -rtuv --delete "$NXSDK" "host-1:NxSDK"
}
pull_host1() {
rsync -rtuv --exclude=*.pyc "host-1:nengo-loihi/" "$LOIHI" rsync -rtuv "host-1:NxSDK" "$NXSDK"
}

These functions are placed in the ~/.bashrc file and executed at a terminal with

push_host1
pull_host1

### Remote editing with SSHFS¶

If you primarily work with a single remote superhost, SSHFS is a good option that allows you to mount a remote filesystem to your local machine, meaning that you manipulate files as you normally would on your local machine, but those files will actually exist on the remote machine. SSHFS ensures that change you make locally are efficiently sent to the remote.

SSHFS is available from most package managers, including apt and brew.

To mount a remote directory to your local machine, create a directory to mount to, then call sshfs to mount it.

mkdir -p <mount point>
sshfs -o allow_other,defer_permissions <host short name>:<remote directory> <mount point>

When you are done using the remote files, unmount the mount point.

fusermount -u <mount point>

Note

If fusermount is not available and you have sudo access, you can also unmount with

sudo umount <mount point>

As with rsync, since you may do these commands frequently, it can save time to make a short bash function. The following example functions mount and unmount the host-2 ~/loihi directory to the local machine’s ~/remote/host-2 directory.

mount_host2() {
mkdir -p ~/remote/host-2
sshfs host-2:loihi ~/remote/host-2
}
unmount_host2() {
fusermount -u ~/remote/host-2
}

## Superhost¶

### Plotting¶

If you are generating plots with Matplotlib on the superhost or host, you may run into issues due to there being no monitor attached to those machines (i.e., they are “headless”). Rather than plotting to a screen, you can instead save plots as files with plt.savefig. You will also need to configure Matplotlib to use a headless backend by default.

The easiest way to do this is with a matplotlibrc file.

mkdir -p ~/.config/matplotlib
echo "backend: Agg" >> ~/.config/matplotlib/matplotlibrc

### IPython / Jupyter¶

If you want to use the IPython interpreter or the Jupyter notebook on a superhost (e.g., the INRC superhost), you may run into issues due to the network file system (NFS), which does not work well with how IPython and Jupyter track command history. You can configure IPython and Jupyter to instead store command history to memory only.

To do this, start by generating the configuration files.

jupyter notebook --generate-config
ipython profile create

Then add a line to three files to configure the command history for NFS.

echo "c.NotebookNotary.db_file = ':memory:'" >> ~/.jupyter/jupyter_notebook_config.py
echo "c.HistoryAccessor.hist_file = ':memory:'" >> ~/.ipython/profile_default/ipython_config.py
echo "c.HistoryAccessor.hist_file = ':memory:'" >> ~/.ipython/profile_default/ipython_kernel_config.py

### Slurm cheatsheet¶

Most Loihi superhosts use Slurm to schedule and distribute jobs to Loihi hosts. Below are the commands that Slurm makes available and what they do.

sinfo

Check the status (availability) of connected hosts.

squeue

Check the status of your jobs.

scancel <jobid>

scancel --user=<username>

sudo scontrol update nodename="<nodename>" state="idle"

Mark a Loihi host as “idle”, which places it in the pool of available hosts to be used. Use this when a Loihi host that was down comes back up.

Note

This should only be done by a system administrator.

### Use Slurm by default¶

Most superhosts use Slurm to run models on the host. Normally you can opt in to executing a command with

SLURM=1 my-command

However, you will usually want to use Slurm, so to switch to an opt-out setup, open your shell configuration file in a text editor (usually ~/.bashrc), and add the following line to the end of the file.

export SLURM=1

Once making this change you can opt out of using Slurm by executing a command with

SLURM=0 my-command

### Running large models¶

Normally you do not need to do anything other than setting the SLURM environment variable to run a model on Slurm. However, in some situation Slurm may kill your job due to long run times or other factors.

Custom Slurm partitions can be used to run your job with different sets of restrictions. Your system administrator will have to set up the partition. You can see a list of all partitions and nodes with sinfo.

To run a job with the loihiinf partition, set the environment variable PARTITION. For example, you can run bigmodel.py using this partition with

PARTITION=loihiinf python bigmodel.py

Similarly, if you wish to use a particular board (called a “node” in Slurm), set the BOARD environment variable. For example, to run model.py on the loihimh board, do

BOARD=loihimh python model.py

format: plain


text:
# Weierstrass Approximation Theorem

## Theorem

Let $f$ be a real function which is continuous on the closed interval $\Bbb I$.

Then $f$ can be uniformly approximated on $\Bbb I$ by a polynomial function to any given degree of accuracy.

## Also known as

This result is also seen referred to as Weierstrass's theorem, but as there are a number of results bearing Karl Weierstrass's name, it makes sense to be more specific.

## Source of Name

This entry was named for Karl Weierstrass.

## Historical Note

The Weierstrass Approximation Theorem has been demonstrated to have far-reaching and important effects in every aspect of the field of analysis.

It has also been given a significant generalisation by Marshall Harvey Stone.


text:
# Clojure Numerics, Part 6 - More Linear Algebra Fun with Least Squares

December 27, 2017

This time, we look at a few important variants of the least squares problem: least squares with equality constraints, and generalized least squares. They show how we can solve problems where there is some norm to be minimized with additional constraints to be satisfied. A very cool and useful thing to know!

Before I continue, a few reminders:

• Include Neanderthal library in your project to be able to use ready-made high-performance linear algebra functions.
• Read articles in the introductory Clojure Linear Algebra Refresher series.
• This is the sixth article in a more advanced series. It is a good idea to read them all in sequence if this is your first contact with this kind of software.

The namespaces we'll use:

(require '[uncomplicate.neanderthal
[native :refer [dge dv]]
[core :refer [mm mv]]
[linalg :refer [lse gls]]])

## Least Squares With Equality Constraints

We have a system $$Ax=c$$ that we need to find the least squares solution, as in the last article. But, additionally, we have the constraint $$Bx=d$$. There is much theory connected to how to do this in a stable way, but the programmer in me just want to skip it and show you the code. So, here it is. In Clojure, just use the lse function!

(let [a (dge 6 4 [-0.57 -1.28 -0.39 0.25
-1.93 1.08 -0.31 -2.14
2.30 0.24 0.40 -0.35
-1.93 0.64 -0.66 0.08
0.15 0.30 0.15 -2.13
-0.02 1.03 -1.43 0.50]
{:layout :row})
b (dge 2 4 [1 0 -1 0
0 1 0 -1]
{:layout :row})
c (dv [-1.50 -2.14 1.23 -0.54 -1.68 0.82])
d (dv [7 2])]
(lse a b c d))

#RealBlockVector[double, n:4, offset: 0, stride:1]
[   2.57    1.12   -4.43   -0.88 ]

Our solution is $$x=(2.57, 1.12, -4.43, -0.88)$$.

Is it correct? Let's see if it satisfies the constraint:

(mv (dge 2 4 [1 0 -1 0
0 1 0 -1]
{:layout :row})
(dv [2.57 1.12 -4.43 -0.88 ]))

#RealBlockVector[double, n:2, offset: 0, stride:1]
[   7.00    2.00 ]

Yes!

I won't bother you longer with this, other than repeat that the search key phrase for this kind of problems is Least Sqares with Equality Constraints (LSE) and that the lse and its impure cousin lse! is what you need in those cases.

## Generalized Least Squares

Here's another stock problem from the Least Squares family. This one is connected to problems commonly found in statistics, data analysis, machine learning, etc. such as regularization, weighting and similar stuff.

Let's see this case similar to the stock problem solvable by least squares: $$Ax+w=c$$. We added $$w$$ to the basic equation. It can describe, for example, normally distributed noise in the system, with zero mean, and covariance matirx $$\sigma^2 W$$, where $$W=BB^T$$. Now, we could have proceeded by doing some matrix algebra and concluding that the solution can be found by minimizing $${\lVert B^{-1}(Ax-c) \rVert}_2$$. We know (since we read that in the previous articles in this series! :) that computing and working with inverse matrices is bad, bad, bad: a slow procedure that often gives wrong answers.

It was shown by smart mathematicians that this problem is equivalent to the problem of minimizing the norm $${\lVert y^Ty \rVert}$$ with the constraint $$Ax+By=c$$. Fair enough, but how does it help us?

That problem is called the Generalized Least Sqares problem (GLS) and it is numerically stable! It is done by a flavor of QR factorization that is too exotic for us to look in detail. The neat thing is that we don't have to, because Neanderthal comes with the function that does everything for us: gls.

(let [a (dge 4 3 [-0.57 -1.28 -0.39
-1.93 1.08 -0.31
2.30 0.24 -0.40
-0.02 1.03 -1.43]
{:layout :row})
b (dge 4 4 [0.5 0 0 0
0 1 0 0
0 0 2 0
0 0 0 5]
{:layout :row})
c (dv [1.32 -4 5.52 3.24])]
(gls a b c))

'(#RealBlockVector(double  n:3  offset: 0  stride:1)
(   1.99   -1.01   -2.99 )
#RealBlockVector(double  n:4  offset: 0  stride:1)
(  -0.00   -0.00   -0.00    0.01 )
)

The solution is given as a Clojure vector; a pair containing $$x$$ and $$y$$. In this example, GLS solution is $$x=(1.99, -1.01, -2.99), y=(0, 0, 0, 0.01)$$ (roughly; here, only the first two decimals are printed).

## Enough of the least squares

Huh. This time I managed to be brief. For better or for worse :)

There are quite a few problems that can be described and solved as one of the variants of least squares. I hope that these typical flavors that I showed can give you enough foundations to find your own way when there's something more exotic to be solved. Don't think that what we'd covered is a small thing, though. Quite a lot of problems could be solved with : ls, lse, and gls. A lot!

Maybe it's time to give some attention to eigen problems now. Remember the Clojure Linear Algebra Refresher (2) - Eigenvalues and Eigenvectors article from the Linear Algebra Refresher series? In the upcoming articles, we'll look into it a bit more seriously.

Clojure Numerics, Part 6 - More Linear Algebra Fun with Least Squares - December 27, 2017 - Dragan Djuric


text:
# How is the set $\{x:F'(x)=\infty\}$ characterized?

I am asking a generalization of the following problem. Given a measure zero set $E$, we can find $f$ so that $F(x)=\int^x_0 f(t) dt$ has $F'(x)=\infty$ on $E$. In particular, by outer regularity, let $V_n \supseteq E$ be open and $|V_n|<2^{-n}$. Then $f(x)=\sum_{n=1} \chi_{V_n}$ does the trick.

But this means that $F'(x)=\infty$ on a little more than $E$, possibly, since we approximated $E$ by a Borel set. Is it possible to find $f$ such that $F'(x)=\infty$ precisely on $E$, given any $E$ with measure zero? Can we characterize the set $\{x:F'(x)=\infty\}$?

-
No, it's not possible. Philosophically speaking, you never get anything non-Borel from a continuous function such as your $F$. But I don't have time to write a detailed proof right now. You may want to read the 2nd reference in the article on Zahorski's theorem, and a more recent paper by Fowler and Preiss. – user31373 Aug 7 '12 at 20:15


text:
A body on the surface of the earth weighs 400 lb. The radius of the earth is about 4,000 miles. If this same body were placed on a 4000 mile high tower, it's weight would be:

a) 100 lbs
b) 200lbs
c) 800 lbs
d) 400 lbs
e) 16000 lbs

### Get this answer with Chegg Study

Practice with similar questions
Q:
A body on the surface of the earth weighs 400 lb. the radius of the earth is about 4,000 miles. If this same body were placed on a 4,000 mile high tower, its weight would bea) 100 lbsb) 200 lbsc) 400 lbsd) 800 lbse) 16,000 lbs


text:
# Recommend caliper ≈ $50 (±) ## Main Question or Discussion Point apologies to all if this post is considered inappropriate: i do not know a mechanical engineer nor have i ever needed to own a caliper. i seek a recommendation for a caliper of any design strategy that costs ≈$50 (±) and is the best-value-for-the-dollar.

Related Mechanical Engineering News on Phys.org
Nidum
Gold Member
Electronic , vernier or dial type and what size range ?

If you don't know for certain then just tell us what you want to use it for .

Last edited:
size range ≤ ½"

i guess dial type; is there an advantage to vernier type? is digital in this price range less accurate or less robust ?

i need to measure cylinder diameters to the nearest 1/64" .. likely less than 15/32"

Nidum
Gold Member
You can buy a mechanical vernier with 1/64 readings but best option would be a low cost general purpose electronic caliper with fractions display option .

Google the actual words ' electronic caliper with fractions ' .

There are very expensive ones and very cheap ones but you should find something of satisfactory quality within your price range .

Nidum
Gold Member
As for size range : 150mm/6inch is a good choice for general purpose use .

You can buy a mechanical vernier with 1/64 readings but best option would be a low cost general purpose electronic caliper with fractions display option .

Google the actual words ' electronic caliper with fractions ' .

There are very expensive ones and very cheap ones but you should find something of satisfactory quality within your price range .
can you identify a particular brand and model because i'm a stranger in a strange land.
would i realize an advantage by using the mechanical vernier ?

what do you think of these three:

http://www.ebay.com/itm/ANTIQUE-STARRETT-USA-OUTSIDE-4-INCH-CALIPER-TOOL-VERY-GOOD-CONDITION-/272211383337?hash=item3f61101429:g:UsYAAOSwu1VW2Kjq

Last edited:
Any auto parts store ought to carry them.. Cheap ones can be had for $15, some are plastic, but still measure to within .005", Metal ones will be accurate to .001", some will have fractional display options to 1/128th or at least 1/64th... the one I have is branded "Storm".. handy thing to have in your pocket, and cheap enough to replace.. If you want a GOOD brand, I would look for a used Mitutoyo,.. they're among the best out there, and I had one I bought used 10 years ago and it's been really good for me Nidum Science Advisor Gold Member One of these is what I had in mind : http://www.arceurotrade.co.uk/Catalogue/Measurement/Calipers/Digital-Vernier-Calipers---150mm-300mm-6-12 [Broken] They are available under different brand names in most countries . Last edited by a moderator: billy_joule Science Advisor If you're left handed I'd recommend getting left handed calipers, it's worth the extra$$for the ease of use. Nidum Science Advisor Gold Member @ billy_joule Yes indeed - though most people think that they don't exist ........ http://www.arceurotrade.co.uk/Catalogue/Measurement/Calipers/Left-Hand-Digital-Vernier-Calipers---Not-just-for-Lefties [Broken] You can get left handed micrometers as well . Last edited by a moderator: billy_joule Science Advisor @ billy_joule Yes indeed - though most people think that they don't exist ........ I'd had RH calipers for almost 10 years before I found out LH ones existed, it was quite a revelation! I hate going back to my RH ones now. Nidum Science Advisor Gold Member Reminds me of all those workshop pranks where the old hands sent raw apprentices to the stores for a left handed screw driver and a long weight ....... Last edited: I'm a southpaw, and I've NEVER seen left handed calipers.. I'm at least somewhat ambidextrous so I do alright with standard calipers.. Now measuring tapes, I'd love to see them in left handed models! Last edited by a moderator: One of these is what I had in mind : http://www.arceurotrade.co.uk/Catalogue/Measurement/Calipers/Digital-Vernier-Calipers---150mm-300mm-6-12 [Broken] They are available under different brand names in most countries . follows is arceurotrade company's response to my request to purchase ... Dear Mr. Weissman, Many thanks for your email. Unfortunately, due to our product liability insurance we are not permitted to ship any electrical based products to the USA. We are more than happy to consider orders for mechanical based products for delivery to the USA. Little Machine Shop may be able to assist you with your purchase requirements. Kind Regards Your Enquiries Team at Arc Euro Trade Ltd. • i'm now back to square-1: what are the brands that i should buy - in the ≈$50 (±) price-range - used or new - in the u.s. that offer precision, durability, and value-for-the-dollar?

• i seek your opinion of this - i mention again that my purpose is to measure cylinder diameters to 1/64" -
http://www.ebay.com/itm/Vintage-Browne-Sharpe-Caliper-in-Original-Box-/351731402722?hash=item51e4d347e2:g:S-MAAOSwInxXM3oS

Last edited by a moderator:
Nidum
Gold Member

I'm a tool and die maker by day, so use my vernier constantly. IMO mitutoyo verniers' have gone down in quality in the past 5-10 years, I have some from the '80's that are in better shape than my new ones. If your dial (mitutoyo) has a little hole in the back of the face to re-set the indicator its a good one! Hang on to it!

Insize is a good up and coming brand it digital and dial, were using a lot more of their stuff.

For the OP, If it were me and you need it just for odds and ends stuff I would be hitting amazon.com and picking up one for 20 bucks with free shipping.
If you get a digital make sure it either has a button to turn it off or an auto off. I have a couple brown and sharpe ones that stay on constantly and battery life sucks.

1800 reviews say its pretty good for the price, at 18bucks its not much of a gamble
https://www.amazon.com/dp/B000GSLKIW/?tag=pfamazon01-20

Last edited by a moderator:


text:
### Archive

Archive for November, 2012

## Extending Commodity time series

November 22, 2012 1 comment

I want to follow up with Extending Gold time series post by showing how we can extend Commodity time series. Most Commodity ETFs began trading in 2006, please see the List of Commodity ETFs page. I will use DBC – PowerShares DB Commodity Fund, one on the most liquid Commodity ETFs as my proxy for Commodity time series.

Unfortunately the historical data for the underlying index for the DBC – PowerShares DB Commodity Fund is not publicly available. Instead I will use Thomson Reuters/Jefferies CRB Index to extend DBC – PowerShares DB Commodity Fund time series.

I created a helper function get.CRB() function in data.r at github to download and extract index information. Please note that I used gdata package to read index data and you might need local Perl available to make it work. There is a great tutorial on setting everything up at How to read an excel file (dot xls and dot xlsx) into a data frame with r: gdata and perl post.

Now let’s plot the CRB index and a few Commodity ETFs side by side to visual check the similarity:

###############################################################################
# Load Systematic Investor Toolbox (SIT)
# https://systematicinvestor.wordpress.com/systematic-investor-toolbox/
###############################################################################
setInternet2(TRUE)
con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
source(con)
close(con)

#*****************************************************************
# Load historical data
#******************************************************************
load.packages('quantmod')
CRB = get.CRB()

tickers = spl('GSG,DBC')
getSymbols(tickers, src = 'yahoo', from = '1970-01-01')

#*****************************************************************
# Compare different indexes
#******************************************************************
out = na.omit(merge(Ad(CRB), Ad(GSG), Ad(DBC)))
colnames(out) = spl('CRB,GSG,DBC')
temp = out / t(repmat(as.vector(out[1,]),1,nrow(out)))

# Plot side by side
layout(1:2, heights=c(4,1))
plota(temp, ylim=range(temp))
plota.lines(temp[,1],col=1)
plota.lines(temp[,2],col=2)
plota.lines(temp[,3],col=3)
plota.legend(colnames(temp),1:3)

# Plot correlation table
temp = compute.cor(temp / mlag(temp)- 1, 'pearson')
temp[] = plota.format(100 * temp, 0, '', '%')
plot.table(temp)

The CRB index performed similar to Commodity ETFs; its not the exact replica, but will do given the lack of public historical data for Commodity index.

Next let’s extend the DBC’s time series with CRB index and use it for a simple equal weight strategy:

	#*****************************************************************
# Create simple equal weight back-test
#******************************************************************
tickers = spl('GLD,DBC,TLT')

data <- new.env()
getSymbols(tickers, src = 'yahoo', from = '1980-01-01', env = data, auto.assign = T)
for(i in ls(data)) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)

# extend Gold and Commodity time series
data$GLD = extend.GLD(data$GLD)
data$DBC = extend.data(data$DBC, get.CRB(), scale=T)

bt.prep(data, align='remove.na')

#*****************************************************************
# Code Strategies
#******************************************************************
prices = data$prices n = ncol(prices) # find period ends period.ends = endpoints(prices, 'months') period.ends = period.ends[period.ends > 0] models = list() #***************************************************************** # Equal Weight #****************************************************************** data$weight[] = NA
data$weight[period.ends,] = ntop(prices[period.ends,], n) models$equal.weight = bt.run.share(data, clean.signal=F)

#*****************************************************************
# Create Report
#******************************************************************
plotbt.custom.report.part1(models)
plotbt.custom.report.part2(models)

By extending DBC’s and GLD’s time series, the back-test now goes back to the July 2002, the TLT’s first trading date.

To view the complete source code for this example, please have a look at the bt.extend.DBC.test() function in bt.test.r at github.

Advertisements
Categories: R

## Regime Detection Pitfalls

Today, I want to address some questions that I was getting regarding the Regime Detection post. In the Regime Detection post I showed an example based on the simulated data, and some of you tried to apply this example to actual stocks.

There is one big problem that you have to be aware in designing the trading system that uses Regime Detection. Let me demonstrate with a simple example.

Let’s say we used 10 years of daily returns of SPY to train our model (2 regimes: Bull and Bear). Next let’s test a strategy that goes long in the Bull regime and stays in cash in the Bear regime using the out of sample data. I thank you Michael Guan at systematicedge.wordpress.com for asking this question.

###############################################################################
# Load Systematic Investor Toolbox (SIT)
# https://systematicinvestor.wordpress.com/systematic-investor-toolbox/
###############################################################################
setInternet2(TRUE)
con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
source(con)
close(con)

#*****************************************************************
# Load historical data
#******************************************************************
load.packages('quantmod')
data <- new.env()
getSymbols('SPY', src = 'yahoo', from = '1970-01-01', env = data, auto.assign = T)
data$SPY = adjustOHLC(data$SPY, use.Adjusted=T)
bt.prep(data)

#*****************************************************************
# Setup
#******************************************************************
nperiods = nrow(data$prices) models = list() rets = ROC(Ad(data$SPY))
rets[1] = 0

# use 10 years: 1993:2002 for training
in.sample.index = '1993::2002'
out.sample.index = '2003::'

in.sample = rets[in.sample.index]
out.sample = rets[out.sample.index]
out.sample.first.date = nrow(in.sample) + 1

#*****************************************************************
# Fit Model
#******************************************************************
load.packages('RHmm')
fit = HMMFit(in.sample, nStates=2)

# find states
states.all = rets * NA
states.all[] = viterbi(fit, rets)$states #***************************************************************** # Code Strategies #****************************************************************** data$weight[] = NA
data$weight[] = iif(states.all == 1, 0, 1) data$weight[in.sample.index] = NA
models$states.all = bt.run.share(data) # create plot plotbt.custom.report.part1(models)  These results are too good to be true for such a simple strategy. You have tripped your investment in 10 years with minimum draw-down. So there is the problem? The Viterbi algorithm uses all data to compute the most likely sequence of states. I.e. In the code above, we passed the returns for all dates to the Viterbi algorithm. Hence, for example, the Viterbi algorithm might use data in 2008 to determine states in 2007. This is a big problem, that we can address using a rolling window approach. Unfortunately, the accuracy of our predictions and strategy returns go down.  #***************************************************************** # Find problem - results are too good #****************************************************************** # The viterbi function need to see all data to compute the most likely sequence of states # http://en.wikipedia.org/wiki/Viterbi_algorithm # We can use expanding window to determine the states states.win1 = states.all * NA for(i in out.sample.first.date:nperiods) { states.win1[i] = last(viterbi(fit, rets[1:i])$states)
if( i %% 100 == 0) cat(i, 'out of', nperiods, '\n')
}

# Or we can refit model over expanding window as suggested in the
# Regime Shifts: Implications for Dynamic Strategies by M. Kritzman, S. Page, D. Turkington
# Out-of-Sample Analysis, page 8
initPoint = fit$HMM states.win2 = states.all * NA for(i in out.sample.first.date:nperiods) { fit2 = HMMFit(rets[2:i], nStates=2, control=list(init='USER', initPoint = initPoint)) initPoint = fit2$HMM
states.win2[i] = last(viterbi(fit2, rets[2:i])$states) if( i %% 100 == 0) cat(i, 'out of', nperiods, '\n') } #***************************************************************** # Plot States #****************************************************************** layout(1:3) col = col.add.alpha('white',210) plota(states.all[out.sample.index], type='s', plotX=F) plota.legend('Implied States based on all data', x='center', bty='o', bg=col, box.col=col,border=col,fill=col,cex=2) plota(states.win1[out.sample.index], type='s') plota.legend('Implied States based on rolling window', x='center', bty='o', bg=col, box.col=col,border=col,fill=col,cex=2) plota(states.win2[out.sample.index], type='s') plota.legend('Implied States based on rolling window(re-fit)', x='center', bty='o', bg=col, box.col=col,border=col,fill=col,cex=2)  The states change more often using the rolling window approach. We probably can add a filter that allows the state change only it is confirmed by certain number of observations. I.e. delayed state change, let’s keep it for another post. Let’s just see the impact of the rolling window approach on the back-test results:  #***************************************************************** # Code Strategies #****************************************************************** data$weight[] = NA
data$weight[] = iif(states.win1 == 1, 0, 1) data$weight[in.sample.index] = NA
models$states.win1 = bt.run.share(data) data$weight[] = NA
data$weight[] = iif(states.win2 == 1, 0, 1) data$weight[in.sample.index] = NA
models$states.win2 = bt.run.share(data) #***************************************************************** # Create report #****************************************************************** plotbt.custom.report.part1(models)  As expected, it is quite hard to apply Regime Detection to stocks. To view the complete source code for this example, please have a look at the bt.regime.detection.pitfalls.test() function in bt.test.r at github. Categories: R ## Simulating Multiple Asset Paths in R I recently came across the Optimal Rebalancing Strategy Using Dynamic Programming for Institutional Portfolios by W. Sun, A. Fan, L. Chen, T. Schouwenaars, M. Albota paper that examines the cost of different rebablancing methods. For example, one might use calendar rebalancing: i.e. rebalance every month / quarter / year. Or one might use threshold rebalancing: i.e. rebalance only if asset weights break out from a certain band around the policy mix, say 3%. To investigate the cost of the different rebalancing methods, authors run 10,000 simulations. Today, I want to show how to simulate asset price paths given the expected returns and covariances. I will assume that prices follow the Geometric Brownian Motion. Also I will show a simple application of Monte Carlo option pricing. In the next post I will evaluate the cost of different rebalancing methods. Let’s assume that a stock price can be described by the stochastic differential equation: $dS_{t}=\mu S_{t}dt+\sigma S_{t}dW_{t}$ where $\mu$ is the expected annual return of the underlying asset, $\sigma$ is the annualized volatility, and $W_{t}$ is a Brownian Motion. This stochastic differential equation has the following solution $S_{t}=S_{0}exp\left(\left(\mu - \frac{\sigma^{2}}{2}\right)t+\sigma W_{t}\right)$ which I implemented in the asset.paths() function. The asset.paths() function is based on the Simulating Multiple Asset Paths in MATLAB code in Matlab. asset.paths <- function(s0, mu, sigma, nsims = 10000, periods = c(0, 1) # time periods at which to simulate prices ) { s0 = as.vector(s0) nsteps = len(periods) dt = c(periods[1], diff(periods)) if( len(s0) == 1 ) { drift = mu - 0.5 * sigma^2 if( nsteps == 1 ) { s0 * exp(drift * dt + sigma * sqrt(dt) * rnorm(nsims)) } else { temp = matrix(exp(drift * dt + sigma * sqrt(dt) * rnorm(nsteps * nsims)), nc=nsims) for(i in 2:nsteps) temp[i,] = temp[i,] * temp[(i-1),] s0 * temp } } else { require(MASS) drift = mu - 0.5 * diag(sigma) n = len(mu) if( nsteps == 1 ) { s0 * exp(drift * dt + sqrt(dt) * t(mvrnorm(nsims, rep(0, n), sigma))) } else { temp = array(exp(as.vector(drift %*% t(dt)) + t(sqrt(dt) * mvrnorm(nsteps * nsims, rep(0, n), sigma))), c(n, nsteps, nsims)) for(i in 2:nsteps) temp[,i,] = temp[,i,] * temp[,(i-1),] s0 * temp } } }  Next let’s visualize some simulation asset paths: ############################################################################### # Load Systematic Investor Toolbox (SIT) # https://systematicinvestor.wordpress.com/systematic-investor-toolbox/ ############################################################################### setInternet2(TRUE) con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb')) source(con) close(con) #***************************************************************** # Plot some price paths #****************************************************************** S = c(100,105) X = 98 Time = 0.5 r = 0.05 sigma = c(0.11,0.16) rho = 0.63 N = 10000 # Single Asset for 10 years periods = 0:10 prices = asset.paths(S[1], r, sigma[1], N, periods = periods) # plot matplot(prices[,1:100], type='l', xlab='Years', ylab='Prices', main='Selected Price Paths') # Multiple Assets for 10 years periods = 0:10 cov.matrix = sigma%*%t(sigma) * matrix(c(1,rho,rho,1),2,2) prices = asset.paths(S, c(r,r), cov.matrix, N, periods = periods) # plot layout(1:2) matplot(prices[1,,1:100], type='l', xlab='Years', ylab='Prices', main='Selected Price Paths for Asset 1') matplot(prices[2,,1:100], type='l', xlab='Years', ylab='Prices', main='Selected Price Paths for Asset 2')  Next, let’s look at examples of Monte Carlo option pricing using asset.paths() function in random.r at github.  #***************************************************************** # Price European Call Option #****************************************************************** load.packages('fOptions') # Black–Scholes GBSOption(TypeFlag = "c", S = S[1], X = X, Time = Time, r = r, b = r, sigma = sigma[1]) # Monte Carlo simulation N = 1000000 prices = asset.paths(S[1], r, sigma[1], N, periods = Time) future.payoff = pmax(0, prices - X) discounted.payoff = future.payoff * exp(-r * Time) # option price mean(discounted.payoff) #***************************************************************** # Price Asian Call Option #****************************************************************** load.packages('fExoticOptions') Time = 1/12 # Approximation TurnbullWakemanAsianApproxOption(TypeFlag = "c", S = S[1], SA = S[1], X = X, Time = Time, time = Time, tau = 0 , r = r, b = r, sigma = sigma[1]) # Monte Carlo simulation N = 100000 periods = seq(0,Time,1/360) n = len(periods) prices = asset.paths(S[1], r, sigma[1], N, periods = periods) future.payoff = pmax(0, colSums(prices)/n - X) discounted.payoff = future.payoff * exp(-r * Time) # option price mean(discounted.payoff) #***************************************************************** # Price Basket Option #****************************************************************** Time = 0.5 # Approximation TwoRiskyAssetsOption(TypeFlag = "cmax", S1 = S[1], S2 = S[2], X = X, Time = Time, r = r, b1 = r, b2 = r, sigma1 = sigma[1], sigma2 = sigma[2], rho = rho) # Monte Carlo simulation N = 100000 cov.matrix = sigma%*%t(sigma) * matrix(c(1,rho,rho,1),2,2) prices = asset.paths(S, c(r,r), sigma = cov.matrix, N, periods = Time) future.payoff = pmax(0, apply(prices,2,max) - X) discounted.payoff = future.payoff * exp(-r * Time) # option price mean(discounted.payoff) #***************************************************************** # Price Asian Basket Option #****************************************************************** Time = 1/12 # Monte Carlo simulation N = 10000 periods = seq(0,Time,1/360) n = len(periods) prices = asset.paths(S, c(r,r), sigma = cov.matrix, N, periods = periods) future.payoff = pmax(0, colSums(apply(prices,c(2,3),max))/n - X) discounted.payoff = future.payoff * exp(-r * Time) # option price mean(discounted.payoff)  Please note that Monte Carlo option pricing requireies many simulations to converge to the option price. It takes longer as we increase number of simulations or number of time periods or number of assets. On the positive side, it provides a viable alternative to simulating difficult problems that might not be solved analytically. In the next post I will look at the cost of different rebalancing methods. To view the complete source code for this example, please have a look at the asset.paths.test() function in random.r at github. Categories: R, Risk Measures ## Regime Detection November 1, 2012 7 comments Regime Detection comes handy when you are trying to decide which strategy to deploy. For example there are periods (regimes) when Trend Following strategies work better and there are periods when Mean Reversion strategies work better. Today I want to show you one way to detect market Regimes. To detect market Regimes, I will fit a Hidden Markov Regime Switching Model on the set of simulated data (i.e. Bull / Bear market environments) I will use the excellent example from the Markov Regime Switching Models in MATLAB post and adapt it to R. The idea behind using the Regime Switching Models to identify market states is that market returns might have been drawn from 2 or more distinct distributions. As a base case, for example, we may suppose that market returns are samples from one normal distribution N(mu, sigma) i.e. Returns = mu + e, e ~ N(0, sigma)  Next we may suppose that market returns are samples from two normal distributions (i.e. returns during Bull market may be ~ N(mu.Bull, sigma.Bull) and returns during Bear market may be N(mu.Bear , sigma.Bear) i.e. Returns = mu + e, e ~ N(0, sigma) mu = mu.Bull for Bull regime and mu.Bear for Bear regime and sigma= sigma.Bull for Bull regime and sigma.Bear for Bear regime  Fortunately we do not have to fit regimes by hand, there is the RHmm package for Hidden Markov Models at CRAN that uses the Baum-Welch algorithm to fit Hidden Markov Models. Next, let follow the steps from the Markov Regime Switching Models in MATLAB post. ############################################################################### # Load Systematic Investor Toolbox (SIT) # https://systematicinvestor.wordpress.com/systematic-investor-toolbox/ ############################################################################### setInternet2(TRUE) con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb')) source(con) close(con) #***************************************************************** # Generate data as in the post #****************************************************************** bull1 = rnorm( 100, 0.10, 0.15 ) bear = rnorm( 100, -0.01, 0.20 ) bull2 = rnorm( 100, 0.10, 0.15 ) true.states = c(rep(1,100),rep(2,100),rep(1,100)) returns = c( bull1, bear, bull2 ) # find regimes load.packages('RHmm') y=returns ResFit = HMMFit(y, nStates=2) VitPath = viterbi(ResFit, y) #Forward-backward procedure, compute probabilities fb = forwardBackward(ResFit, y) # Plot probabilities and implied states layout(1:2) plot(VitPath$states, type='s', main='Implied States', xlab='', ylab='State')

matplot(fb$Gamma, type='l', main='Smoothed Probabilities', ylab='Probability') legend(x='topright', c('State1','State2'), fill=1:2, bty='n')  The first chart shows states (1/2) determined by the model. The second chart shows the probability of being in each state. Next, let’s generate some additional data and see if the model is able to identify the regimes  #***************************************************************** # Add some data and see if the model is able to identify the regimes #****************************************************************** bear2 = rnorm( 100, -0.01, 0.20 ) bull3 = rnorm( 100, 0.10, 0.10 ) bear3 = rnorm( 100, -0.01, 0.25 ) y = c( bull1, bear, bull2, bear2, bull3, bear3 ) VitPath = viterbi(ResFit, y)$states

#*****************************************************************
# Plot regimes
#******************************************************************
load.packages('quantmod')
data = xts(y, as.Date(1:len(y)))

layout(1:3)
plota.control$col.x.highlight = col.add.alpha(true.states+1, 150) plota(data, type='h', plotX=F, x.highlight=T) plota.legend('Returns + True Regimes') plota(cumprod(1+data/100), type='l', plotX=F, x.highlight=T) plota.legend('Equity + True Regimes') plota.control$col.x.highlight = col.add.alpha(VitPath+1, 150)
plota(data, type='h', x.highlight=T)
plota.legend('Returns + Detected Regimes')

The first 300 observations were used to calibrate this model, the next 300 observations were used to see how the model can describe the new infromation. This model does relatively well in our toy example.

To view the complete source code for this example, please have a look at the bt.regime.detection.test() function in bt.test.r at github.

Categories: R


text:
# For which integers $a,b,c,d$ does $\frac{a}{b}+\frac{c}{d} = \frac{a+c}{b+d}$?

For which $a,b,c,d \in \mathbb{Z}$ does $\frac{a}{b}+\frac{c}{d} = \frac{a+c}{b+d}$?

This is actually the question I meant to ask in a previous question that I asked here.

What about $a,b,c,d \in \mathbb{R}?$

$$\frac{a}{b} + \frac{c}{d} = \frac{a+c}{b+d}$$

We will keep careful track of which quantities may be zero. To start, we know that $b, d, b+d \neq 0$ because they are in the denominator.

Let's multiply both sides by $b\cdot d \cdot (b+d)$:

\begin{align*}ad(b+d) + bc(b+d) &= (a + c)bd\\ abd + ad^2 + b^2c + bcd &= abd + bcd\\ abd + ad^2 + b^2c + bcd - abd - bcd &= 0\\ \fbox{ad^2 + b^2 c = 0} \end{align*}

Here's an example assignment that satisfies this equation: \begin{align*}a &= 5\\b&=3\\c&=-20\\d&=6\\\end{align*}

We can parameterize this set of solutions with three parameters, $\alpha, \beta_1, \beta_2 \in \mathbb{R}$ and three minimal constraints $\beta_1 \neq 0, \beta_2 \neq 0$, and $\beta_1 \neq \beta_2$. For any such triple, a solution is given by \begin{align*} a &\equiv \alpha\\ b &\equiv \frac{1}{\beta_1}\\ c &\equiv -\alpha\frac{\beta_1^2}{\beta_2^2}\\ d &\equiv -\frac{1}{\beta_2}\\ \end{align*}

Conversely, given any solution $\langle a, b, c, d\rangle$ with $\frac{a}{b} + \frac{c}{d} = \frac{a+c}{b+d}$, we can solve uniquely for parameters $\langle \alpha, \beta_1, \beta_2\rangle$ meeting these constraints.

For fun, here is a nomogram that allows you to find such collections of numbers:

You can draw any line passing through, say, the B and C axes, and any line passing through the A and D axes. If the two lines meet at the same point on the central vertical axis, then the corresponding points $(a, b, c, d)$ satisfy the equation.

In the picture, the dashed lines show an example. Here, we see that the equation is satisfied when $a = 4$, $b=3$, $c=-1$, and $d=1.5$.

Also, if you specify any four of the following five variables—a value for $a$, a value for $b$, a value for $c$, a value for $d$, a spot on the central axis— you can "solve" for the remaining point.

• A parametrization of all the solutions would be nice however :-p. – YoTengoUnLCD Jul 17 '16 at 5:55
• How did you draw the nomogram? Did you use pynomo.org ? – Rodrigo de Azevedo Jul 17 '16 at 12:56
• @RodrigodeAzevedo Yes! I used Pynomo. – user326210 Jul 17 '16 at 17:44
• Looks like you also fell under Ron Doerfler's spell ;-) – Rodrigo de Azevedo Jul 17 '16 at 17:56

If you clear the fractions you get $abd+ad^2+b^2c+bcd=abd+bcd$ or $$ad^2+b^2c=0$$ As long as none of the denominators are zero, you are set.

Multiply both sides by $bd(b+d)$ and you get:

\begin{align}ad(d+b) + cb (d+b) = bd(a+c)&\iff\\ ad^2 + abd + cbd + cb^2 = b d(a+c)&\iff\\ ad^2 +cb^2 = 0&\iff\\ ad^2=-cb^2 \end{align}

Given $b,c$, you need $d^2$ to be a factor of $cb^2$ and then $a=-\frac{cb^2}{d^2}$.

Given non-zero integers $m,n,p,q$, with you have:

$$(a,b,c,d)=(pm^2,mq,-pn^2,nq)$$

This gives all solutions, excluding the case where $m=-n$.

The only solution with $\gcd(a,b)=1$ and $\gcd(c,d)=1$ are the cases where $m,n=\pm 1$, and then you are dealing with the solutions:

$$\frac{a}{b}+\frac{-a}{b}=\frac{a+(-a)}{b+b}$$


text:
# trustedcoin

A lightningd plugin that replaces a bitcoind node with trusted public explorers. Great so one can run a c-lightning node on a VPS without having to bother with an expensive bitcoind.

It uses multiple explorers and alternates between them, parses blocks and compares them with block hashes taken from other explorers so if someone is trying to trick you at least your node will crash instead of assuming all is fine.

It is used on @lntxbot and Etleneum.


text:
# Function that equals 1 at x=0, but 0 everywhere else?

1. Mar 10, 2014

### Kepler_

Is there a function that equals 1 at x=0 and equals 0 when x isn't 0 without using a piecewise function? I've been experimenting with limits and derivatives but haven't made much progress.

The closest thing I was able to think of is 1-x/x, which is indeterminate at x=0, and 0 everywhere else.

Thanks!

2. Mar 10, 2014

### jfizzix

How about the limit as $a$ approaches $\infty$ of $f(x) = e^{-a x^{2}}$?

3. Mar 10, 2014

### Kepler_

That'll do it :P
It seems like any real constant greater than 1 works in place of e. Is there a reason why you chose e?

Last edited: Mar 10, 2014
4. Mar 10, 2014

### DrewD

Out of infinitely many possible bases to choose from, $e$ is the only one that makes sense most of the time. In this case it doesn't matter, but why would you choose something else? If you want, you could rewrite this as
$\lim_{a\rightarrow\infty}\left(\frac{1}{e^a}\right)^{x^2}$in which case it is the same as
$\lim_{b\rightarrow0}\left( b\right)^{x^2}$
What matters is the $x^2$ exponent and the continuous, decreasing function as a base. Most mathematicians are used to $e$ being part of that base.

Also, this is sometimes thought of as a limit of normal distributions. These are defined using $e$.

5. Mar 11, 2014

### Number Nine

Is there a reason for this requirement? "Piecewise" is just a notation for describing a function; it isn't actually a "kind" of function. Even your description of your function is piecewise: 1 at zero, and 0 everywhere else.

6. Mar 11, 2014

### PeroK

Two functions are equal if they take the same value at every point. So, regardless of whether you describe a function piecewise or as the limit of a sequence of functions, it's the same function in the end.

7. Mar 11, 2014

### johnqwertyful

I think this demonstrates a fundamental misunderstanding of functions. Functions are a rule associating two sets. That's it. Really, "piecewise function" doesn't really make sense. It's just a way of writing things.

8. Mar 11, 2014

### DrewD

But sometimes piecewise makes calculus hard. I think this is a good point to make though.

9. Mar 12, 2014

### johnqwertyful

But "Piecewise" is only a way of describing a function. There is no such thing as a "piecewise function". f(x)=x if x>0, -x if x≤0. Or f(x)=|x|. Both are the same function. One is written "piecewise" the other isn't.

10. Mar 12, 2014

### DrewD

Agreed.

11. Mar 12, 2014

### paisiello2

Here's another function that would satisfy the requirement but no limits involved:

f(x) = 0^x

12. Mar 12, 2014

### micromass

Staff Emeritus
How would that be defined for negative real numbers?

13. Mar 12, 2014

### paisiello2

Ooops, you're right, I forgot to modify it:

f(x) = 0^(x^2)

14. Mar 12, 2014

### DrewD

Everytime I have encountered it, $0^0$ has been considered an indeterminate form and not 1.

15. Mar 12, 2014

### micromass

Staff Emeritus
I agree, but let's not open that can of worms here. Paisiello2 clearly meant $0^0=1$, so let's keep it at that, even if it's nonstandard.

mod note: not a $0^0$ argument please. All posts on this will be deleted.

Last edited: Mar 13, 2014


text:
# Microeconomics: Theory Through Applications, v. 1.0

by Russell Cooper and A. Andrew John

## 17.5 Discounted Present Value

Discounted present value is a technique used to add dollar amounts over time. We need this technique because a dollar today has a different value from a dollar in the future.

The discounted present value this year of $1.00 that you will receive next year is If the nominal interest rate is 10 percent, then the nominal interest factor is 1.1, so$1 next year is worth $1/1.1 =$0.91 this year. As the interest rate increases, the discounted present value decreases.

More generally, we can compute the value of an asset this year from this formula:

The flow benefit depends on the asset. For a bond, the flow benefit is a coupon payment. For a stock, the flow benefit is a dividend payment. For a fruit tree, the flow benefit is the yield of a crop.

If an asset (such as a bond) yields a payment next year of $10 and has a price next year of$90, then the “flow benefit from asset + price of the asset next year” is $100. The value of the asset this year is then . If the nominal interest rate is 20 percent, then the value of the asset is$100/1.2 = 83.33.

We discount nominal flows using a nominal interest factor. We discount real flows (that is, flows already corrected for inflation) using a real interest factor, which is equal to (1 + real interest rate).

## Key Insights

• If the interest rate is positive, then the discounted present value is less than the direct sum of flows.
• If the interest rate increases, the discounted present value will decrease.

## More Formally

Denote the dividend on an asset in period t as Dt. Define Rt as the cumulative effect of interest rates up to period t. For example, R2 = (1 + r1)(1 + r2). Then the value of an asset that yields Dt dollars in year t is given by

If the interest rate is constant (equal to r), then the one period interest factor is R = 1 + r, and Rt = Rt.

The discounted present value tool is illustrated in Table 17.1 "Discounted Present Value with Different Interest Rates". The number of years (T) is set equal to 5. The table gives the value of the dividends in each year and computes the discounted present values for two different interest rates. For this example, the annual interest rates are constant over time.

Table 17.1 Discounted Present Value with Different Interest Rates

Year Dividend ($) Discounted Present Value with R = 1.05 ($) Discounted Present Value with R = 1.10 (\$)
1 100 100 100
2 100 95.24 90.91
3 90 81.63 74.38
4 120 103.66 90.16
5 400 329.08 273.20
Discounted Present Value 709.61 628.65


text:
# VR Play postion Sync over the network

I am developing a multiplayer game where a normal character controller user and VR user can participate. So there are two modes of the games through a user can join:

1. Normal
2. VR Mode

I am able to manage these modes and user can participate with any mode but the problem is Normal Player(charcter controller) is unable to view VR controller(HTC vive Controllers) as these objects are becomes off on the side of character controller. So how do I view VR controller at normal player that what VR player is doing.

I make an object and tried to imitiate it position and rotation with respect to VR controller but it is not working as Imitating source(VR controllers) are off (at the normal player side but working in VR mode) and their postion is not changing. How do you manage to make multiplayer game with VR + normal player? or what else the solution.

I have managed to solve this problem using these steps after trying hard three four days. These step are given below so that future user don't get stuck in this problem as there are no other comprehensive guide available:

1. Make one Gameobject for Head (simple cube) with Network Identity and network transform

2. Make one Gameobject for(right) Controller (simple cube) with Network Identity and network transform

3. Make one Gameobject for(left) other Controller (simple cube) with network Identity and network transform

4. Make prefab of all above gameobjects.

5. Add all three prefabs into Network Manager (registered Spawnable prefabs list)

6. delete three prefab from scene

7. Added below script(See comments for details) into my vr player and assign respective prefabs and gameobjects

public class VRPlayerCtrl : NetworkTransform {

//source gameobjects head, left and right controller object of htc vive prefab

public GameObject rightContSource;

public GameObject leftContSource;

//prefabs to assign head, left controller, and right controller
public GameObject vrLeftCtrlPrefab;
public GameObject vrRightCtrlPrefab;

GameObject vrLeftCtrl;
GameObject vrRightCtrl;

void Start()
{

Debug.Log("Start of the vr player");

if (isLocalPlayer)
{
//instantiate prefabs
//disabled conroller meshes at VR player side so it cannont be view by him
vrLeftCtrl.GetComponent<MeshRenderer>().enabled = false;
vrRightCtrl.GetComponent<MeshRenderer>().enabled = false;
}
}

//Instantiate on start head and vr controller object so that it can be view by normal players
{
Debug.Log("instantiateing the controller and head object");
vrLeftCtrl = (GameObject)Instantiate(vrLeftCtrlPrefab);
vrRightCtrl = (GameObject)Instantiate(vrRightCtrlPrefab);

// spawn the bullet on the clients
NetworkServer.Spawn(vrLeftCtrl);
NetworkServer.Spawn(vrRightCtrl);
}

void Update()
{
if (!isLocalPlayer)
{
return;
}

//sync pos on network
CmdControllerPositionSync();
}
//sync position on VR controller objects so that VR player movemnts/action can be viewd by normal user
[Command]
public void CmdControllerPositionSync()
{


text:
# --to find ne income percentage

1 answer below »
need to know formula to find net income percentage answer line looks like the following ;
net income percentage = ? = $= % ------ ----------- ?$

what goe where the ? marks are just need to know the categories are . thanks for any help you can give

## 1 Approved Answer

Mark B
3 Ratings, (9 Votes)
The net Income percentage formula is =...

## Recent Questions in Financial Accounting

Submit Your Questions Here !
Copy and paste your question here...
Attach Files


text:
Mathematics education uncovered and recovered

exam papers

Curious case of an error in an exam paper

Here is a question from A level Edexcel Core 1 January 2013 paper: The equation $$(k+3) x^2+6x+k=5$$   where $$k$$ is a constant, has two distinct real solutions for $$x$$. Show that k satisfies $$k^2-2k-24<0$$   Hence find the set of possible values of $$k$$. This is a question […]

Tags: ,