Volume 29, Issue 3

Low-rank matrix approximation is extremely useful in the analysis of data that arises in scientific computing, engineering applications, and data science. However, as data sizes grow, traditional low-rank matrix approximation methods, such as singular value decomposition (SVD) and column pivoting QR decomposition (CPQR), are either prohibitively expensive or cannot provide sufficiently accurate results. A solution is to use randomized low-rank matrix approximation methods such as randomized SVD, and randomized LU decomposition on extremely large data sets. In this paper, we focus on the randomized LU decomposition method. Then we propose a novel randomized LU algorithm, called SubspaceLU, for the fixed low-rank approximation problem. SubspaceLU is based on the sketch of the co-range of input matrices and allows for an arbitrary number of passes of the input matrix, v ≥ 2 {v\geq 2} . Numerical experiments on CPU show that our proposed SubspaceLU is generally faster than the existing randomized LU decomposition, while remaining accurate. Experiments on GPU shows that our proposed SubspaceLU can gain more speedup than the existing randomized LU decomposition. We also propose a version of SubspaceLU, called SubspaceLU_FP, for the fixed precision low-rank matrix approximation problem. SubspaceLU_FP is a post-processing step based on an efficient blocked adaptive rank determination Algorithm 5 proposed in this paper. We present numerical experiments that show that SubspaceLU_FP can achieve close results to SVD but faster in speed. We finally propose a single-pass algorithm based on LU factorization. Tests show that the accuracy of our single-pass algorithm is comparable with the existing single-pass algorithms.

We study the convergence of Langevin-simulated annealing type algorithms with multiplicative noise, i.e. for V : R d → R V\colon\mathbb{R}^{d}\to\mathbb{R} a potential function to minimize, we consider the stochastic differential equation d ⁢ Y t = − σ ⁢ σ ⊤ ⁢ ∇ V ⁢ ( Y t ) ⁢ d ⁢ t + a ⁢ ( t ) ⁢ σ ⁢ ( Y t ) ⁢ d ⁢ W t + a ⁢ ( t ) 2 ⁢ Υ ⁢ ( Y t ) ⁢ d ⁢ t dY_{t}=-\sigma\sigma^{\top}\nabla V(Y_{t})\,dt+a(t)\sigma(Y_{t})\,dW_{t}+a(t)^{2}\Upsilon(Y_{t})\,dt , where ( W t ) (W_{t}) is a Brownian motion, σ : R d → M d ⁢ ( R ) \sigma\colon\mathbb{R}^{d}\to\mathcal{M}_{d}(\mathbb{R}) is an adaptive (multiplicative) noise, a : R + → R + a\colon\mathbb{R}^{+}\to\mathbb{R}^{+} is a function decreasing to 0 and where Υ is a correction term. Allowing 𝜎 to depend on the position brings faster convergence in comparison with the classical Langevin equation d ⁢ Y t = − ∇ V ⁢ ( Y t ) ⁢ d ⁢ t + σ ⁢ d ⁢ W t dY_{t}=-\nabla V(Y_{t})\,dt+\sigma\,dW_{t} . In a previous paper, we established the convergence in L 1 L^{1} -Wasserstein distance of Y t Y_{t} and of its associated Euler scheme Y ¯ t \bar{Y}_{t} to argmin ⁡ ( V ) \operatorname{argmin}(V) with the classical schedule a ⁢ ( t ) = A ⁢ log − 1 / 2 ⁡ ( t ) a(t)=A\log^{-1/2}(t) . In the present paper, we prove the convergence in total variation distance. The total variation case appears more demanding to deal with and requires regularization lemmas.

The paper gives new results for the Milstein scheme of stochastic differential equations. We show that (i) the Milstein scheme holds as a weak approximation in total variation sense and is given by second-order polynomials of Brownian motion without using iterated integrals under non-commutative vector fields; (ii) the accuracy of the Milstein scheme is better than that of the Euler–Maruyama scheme in an asymptotic sense. In particular, we prove d TV ⁢ ( X T ε , X ¯ T ε , Mil , ( n ) ) ≤ C ⁢ ε 2 / n d_{\mathrm{TV}}(X_{T}^{\varepsilon},\bar{X}_{T}^{\varepsilon,\mathrm{Mil},(n)})\leq C\varepsilon^{2}/n and d TV ⁢ ( X T ε , X ¯ T ε , EM , ( n ) ) ≤ C ⁢ ε / n d_{\mathrm{TV}}(X_{T}^{\varepsilon},\bar{X}_{T}^{\varepsilon,\mathrm{EM},(n)})\leq C\varepsilon/n , where d TV d_{\mathrm{TV}} is the total variation distance, X ε X^{\varepsilon} is a solution of a stochastic differential equation with a small parameter 𝜀, and X ¯ ε , Mil , ( n ) \bar{X}^{\varepsilon,\mathrm{Mil},(n)} and X ¯ ε , EM , ( n ) \bar{X}^{\varepsilon,\mathrm{EM},(n)} are the Milstein scheme without iterated integrals and the Euler–Maruyama scheme, respectively. In computational aspect, the scheme is useful to estimate probability distribution functions by a simple simulation without Lévy area computation. Numerical examples demonstrate the validity of the method.

It is known that the block-based version of the bootstrap method can be used for distributional parameter estimation of dependent data. One of the advantages of this method is that it improves mean square errors. The paper makes two contributions. First, we consider the moving blocking bootstrap method for estimation of parameters of the autoregressive model. For each block, the parameters are estimated based on the modified maximum likelihood method. Second, we provide a method for model selection, Vuong’s test and tracking interval, i.e. for selecting the optimal model for the innovation’s distribution. Our analysis provides analytic results on the asymptotic distribution of the bootstrap estimators and also computational results via simulations. Some properties of the moving blocking bootstrap method are investigated through Monte Carlo simulation. This simulation study shows that, sometimes, Vuong’s test based on the modified maximum likelihood method is not able to distinguish between the two models; Vuong’s test based on the moving blocking bootstrap selects one of the competing models as optimal model. We have studied real data, the S&P500 data, and select optimal model for this data based on the theoretical results.

The partitioning algorithm is an iterative procedure that computes explicitly the steady-state probability of a finite Markov chain 𝑋. In this paper, we propose to adapt this algorithm to the case where the state space E := C ∪ D E:=C\cup D is composed of a continuous part 𝐶 and a finite part 𝐷 such that C ∩ D = ∅ C\cap D=\emptyset . In this case, the steady-state probability 𝜋 of 𝑋 is a convex combination of two steady-state probabilities π C \pi_{C} and π D \pi_{D} of two Markov chains on 𝐶 and 𝐷 respectively. The obtained algorithm allows to compute explicitly π D \pi_{D} . If π C \pi_{C} cannot be computed explicitly, our algorithm approximates it by numerical resolution of successive integral equations. Some numerical examples are studied to show the usefulness and proper functioning of our proposal.