Ivan Dimov
email: ivimov@bas.bg

Sylvain Maire
email: sylvain.maire@univ-tln.fr

Jean Michel Sellier
email: jeanmichel.sellier@parallel.bas.bg

A New Walk on Equations Monte Carlo Method for Linear Algebraic Problems ⋆

Keywords: Monte Carlo algorithms, Markov chain, eigenvalue problem, robust algorithms. 1991 MSC: 65C05, 65U05, 65F10, 65Y20

A new Walk on Equations (WE) Monte Carlo algorithm for Linear Algebra (LA) problem is proposed and studied. This algorithm relies on a non-discounted sum of an absorbed random walk. It can be applied for either real or complex matrices. Several techniques like simultaneous scoring or the sequential Monte Carlo method are applied to improve the basic algorithm. Numerical tests are performed on examples with matrices of different size and on systems coming from various applications. Comparisons with standard deterministic or Monte Carlo algorithms are also done.

Introduction

In this paper we analyze applicability and efficiency of the proposed WE algorithm. The algorithm is based on Markov chain Monte Carlo (MC) approach. We focus our consideration to real symmetric matrices, but a wider class of matrices can also be treated.

Many scientific and engineering applications are based on the problems of solving systems of linear algebraic equations. For some applications it is also important to compute directly the inner product of a given vector and the solution vector of a linear algebraic system. For many LA problems it is important to be able to find a preconditioner with relatively small computational complexity. That is why, it is important to have relatively cheap algorithm for matrix inversion. The computation time for very large problems, or for finding solutions in real-time, can be prohibitive and this prevents the use of many established algorithms. Monte Carlo algorithms give statistical estimates of the required solution, by performing random sampling of a random variable, whose mathematical expectation is the desired solution [START_REF] Curtiss | Monte Carlo methods for the iteration of linear operators[END_REF][START_REF] Curtiss | A Theoretical Comparison of the Efficiencies of two classical methods and a Monte Carlo method for Computing one component of the solution of a set of Linear Algebraic Equations[END_REF][START_REF] Dimov | Monte Carlo Methods for Applied Scientists[END_REF][START_REF] Hammersley | Monte Carlo methods[END_REF][START_REF] Sobol | Monte Carlo numerical methods[END_REF][START_REF] Spanier | Monte Carlo Principles and Neutron Transport Problem[END_REF][START_REF] Westlake | A Handbook of Numerical matrix Inversion and Solution of Linear Equations[END_REF]. Let J be the exact solution of the problem under consideration. Suppose it is proved that there exists a random variable θ, such that E{θ} = J. If one can produce N values of θ, i.e. θ 1 , . . . , θ N , then the value θN = ∑ N i=1 θ N can be considered as a Monte Carlo (MC) approximation of J. The probability error of the MC algorithm is defined as the least possible real number R N , for which: P = P r{| ∑ N i=1 θ N -J| ≤ R N }, where 0 < P < 1.

Several authors have presented works on MC algorithms for LA problems and on the estimation of computational complexity [START_REF] Alexandrov | Parallel Quasi-Monte Carlo Methods for Linear Algebra Problems[END_REF][START_REF] Dimov | Monte Carlo Algorithms for Linear Problems[END_REF][START_REF] Dimov | Monte Carlo algorithms: performance analysis for some computer architectures[END_REF][START_REF] Dimov | Parallel resolvent Monte Carlo algorithms for linear algebra problems[END_REF][START_REF] Dimov | Iterative Monte Carlo algorithms for linear algebra problems[END_REF][START_REF] Dimov | A New Highly Convergent Monte Carlo Method for Matrix Computations[END_REF]. There are MC algorithms for

• computing components of the solution vector;

• evaluating linear functionals of the solution;

• matrix inversion, and • computing the extremal eigenvalues.

An overview of all these algorithms is given in [START_REF] Dimov | Monte Carlo Methods for Applied Scientists[END_REF]. The well-known Power method [START_REF] Golub | Matrix computations[END_REF] gives an estimate for the dominant eigenvalue λ 1 . This estimate uses the so-called Rayleigh quotient. In [START_REF] Dimov | Robustness and applicability of Markov chain Monte Carlo algorithms for eigenvalue problems[END_REF] the authors consider bilinear forms of matrix powers, which is used to formulate a solution for the eigenvalue problem. A Monte Carlo approach for computing extremal eigenvalues of real symmetric matrices as a special case of Markov chain stochastic method for computing bilinear forms of matrix polynomials is considered. In [START_REF] Dimov | Robustness and applicability of Markov chain Monte Carlo algorithms for eigenvalue problems[END_REF] the robustness and applicability of the Almost Optimal Monte Carlo algorithm for solving a class of linear algebra problems based on bilinear form of matrix powers (v, A k ). It is shown how one has to choose the acceleration parameter q in case of using Resolvent Power MC. The systematic error is analyzed and it is shown that the convergence can not be better than O

( 1+|q|λn 1+|q|λ n-1
) m , where λ n is the smallest by modulo eigenvalue, q is the characteristic parameter, and m is the number of MC iterations by the resolvent matrix. To construct an algorithm for evaluating the minimal by modulo eigenvalue λ n , one has to consider the following matrix polynomial p k (A) = ∑ ∞ k=0 q k C k m+k-1 A k , where C k m+k-1 are binomial coefficients, and the characteristic parameter q is used as acceleration parameter of the algorithm [START_REF] Dimov | Parallel resolvent Monte Carlo algorithms for linear algebra problems[END_REF][START_REF] Dimov | Parallel computations of eigenvalues based on a Monte Carlo approach[END_REF][START_REF] Dimov | A Power Method with Monte Carlo Iterations[END_REF]. This approach is a discrete analogues of the resolvent analytical continuation method used in the functional analysis [START_REF] Kantorovich | Approximate Methods of Higher Analysis[END_REF]. There are cases when the polynomial becomes the resolvent matrix [START_REF] Dimov | A New Highly Convergent Monte Carlo Method for Matrix Computations[END_REF][START_REF] Dimov | Parallel computations of eigenvalues based on a Monte Carlo approach[END_REF][START_REF] Dimov | Monte Carlo Methods for Applied Scientists[END_REF]. It should be mentioned that the use of acceleration parameter based on the resolvent presentation is one way to decrease the computational complexity. Another way is to apply a variance reduction technique [START_REF] Densmore | Variational variance reduction for particle transport eigenvalue calculations using Monte Carlo adjoint simulation[END_REF] in order to get the required approximation of the solution with a smaller number of operations. The variance reduction technique for particle transport eigenvalue calculations proposed in [START_REF] Densmore | Variational variance reduction for particle transport eigenvalue calculations using Monte Carlo adjoint simulation[END_REF] uses Monte Carlo estimates of the forward and adjoint fluxes. In [START_REF] Rosca | A new Monte Carlo estimator for systems of linear equations[END_REF] an unbiased estimation of the solution of the system of linear algebraic equations is presented. The proposed estimator can be used to find one component of the solution. Some results concerning the quality and the properties of this estimator are presented. Using this estimator the author gives error bounds and constructs confidence intervals for the components of the solution. In [START_REF] Forsythe | Matrix Inversion by a Monte Carlo Method[END_REF] a Monte Carlo algorithm for matrix inversion is proposed and studied. The algorithm is based on the solution of simultaneous linear equations. In our further consideration we will use some results from [START_REF] Forsythe | Matrix Inversion by a Monte Carlo Method[END_REF] and [START_REF] Dimov | Monte Carlo Methods for Applied Scientists[END_REF] to show how the proposed algorithm can be used for approximation of the inverse of a matrix.

Formulation of the Problem: Solving Linear Systems and Matrix Inversion

By A and B we denote matrices of size n × n, i.e., A, B ∈ IR n×n . We use the following presentation of matrices:

A = {a ij } n i,j=1 = (a 1 , . . . , a i , . . . , a n ) t ,
where a i = (a i1 , . . . , a in ), i = 1, . . . , n and the symbol t means transposition.

The following norms of vectors (l 1 -norm):

∥ b ∥=∥ b ∥ 1 = n ∑ i=1 |b i |, ∥ a i ∥=∥ a i ∥ 1 = n ∑ j=1 |a ij |
and matrices

∥ A ∥=∥ A ∥ 1 = max j n ∑ i=1 |a ij | are used.
Consider a matrix A ∈ IR n×n and a vector b = (b 1 , . . . , b n ) t ∈ IR n×1 . The matrix A can be considered as a linear operator A[IR n → IR n ] [START_REF] Kantorovich | Approximate Methods of Higher Analysis[END_REF][START_REF] Dimov | Monte Carlo Methods for Applied Scientists[END_REF], so that the linear transformation

Ab ∈ IR n×1 (1)
defines a new vector in IR n×1 .

Since iterative Monte Carlo algorithms using the transformation (1) will be considered, the linear transformation ( 1) is called iteration. The algebraic transformation (1) plays a fundamental role in the iterative Monte Carlo algorithms [START_REF] Dimov | Monte Carlo Methods for Applied Scientists[END_REF]. Now consider the following three LA problems Pi (i=1,2,3):

Problem P1. Evaluating the inner product

J(u) = (v, x) = n ∑ i=1 v i x i
of the solution x ∈ IR n×1 of the linear algebraic system

Bx = f, where B = {b ij } n i,j=1 ∈ IR n×n is a given matrix; f = (f 1 , . . . , f n ) t ∈ IR n×1 and v = (v 1 , . . . , v n ) t ∈ IR n×1 are given vectors.
It is possible to choose a non-singular matrix M ∈ IR n×n such that M B = I -A, where I ∈ IR n×n is the identity matrix and

M f = b, b ∈ IR n×1 . Then x = Ax + b. (2)
It will be assumed that

(i)     
1. The matrices M and A are both non-singular;

2. |λ(A)| < 1 for all eigenvalues λ(A) of A,
that is, all values λ(A) for which

Ax = λ(A)x
is satisfied. If the conditions (i) are fulfilled, then a stationary linear iterative algorithm [START_REF] Dimov | Monte Carlo Methods for Applied Scientists[END_REF][START_REF] Dimov | Monte Carlo Algorithms for Linear Problems[END_REF] can be used:

x k = Ax k-1 + b, k = 1, 2, . . . (3)
and the solution x can be presented in a form of a von Neumann series

x = ∞ ∑ k=0 A k b = b + Ab + A 2 b + A 3 b + . . . (4)
The stationary linear iterative Monte Carlo algorithm [START_REF] Dimov | Monte Carlo Methods for Applied Scientists[END_REF][START_REF] Dimov | Monte Carlo Algorithms for Linear Problems[END_REF] is based on the above presentation and will be given later on. As a result, the convergence of the Monte Carlo algorithm depends on the truncation error of the series (3) (see, [START_REF] Dimov | Monte Carlo Methods for Applied Scientists[END_REF]).

Problem P2. Evaluating all components x i , i = 1, . . . n of the solution vector.

In this case we deal with the equation ( 2) assuming that A = I -M B, as well as, b = M f . Problem P3. Inverting of matrices, i.e. evaluating of matrix

G = B -1 ,
where B ∈ IR n×n is a given real matrix.

Assumed that the following conditions are fulfilled:

(ii)

     1. The matrix B is non-singular; 2. ||λ(B)| -1| < 1 for all eigenvalues λ(B) of B.
Obviously, if the condition (i) is fulfilled, the solution of the problem P1 can be obtained using the iterations (3).

For problem P3 the following iterative matrix:

A = I -B can be constructed.
Since it is assumed that the conditions (ii) are fulfilled, the inverse matrix G = B -1 can be presented as

G = ∞ ∑ i=0 A i .
For problems P i (i = 1, 2, 3) one can create a stochastic process using the matrix A and vectors b, and possibly v.

Consider an initial density vector p = {p i } n i=1 ∈ IR n , such that p i ≥ 0, i = 1, . . . , n and p αs > 0 when v αs ̸ = 0 p αs = 0 when v αs = 0.

∑ n i=1 p i = 1.
(

) 5
Similarly, the transition density matrix

P = {p ij } n i,j=1 is called permissible to the matrix A = {a ij } n i,j=1 , i.e. P ∈ P A , if     
p α s-1 ,αs > 0 when a α s-1 ,αs ̸ = 0 p α s-1 ,αs = 0 when a α s-1 ,αs = 0.

In this work we will be dealing with with permissible densities. It is obvious that in such a way the random trajectories constructed to solve the problems under consideration never visit z ero elements of the matrix. Such an approach decreases the computational complexity of the algorithms. It is also very convenient when large sparse matrices are to be treated.

We shall use the so-called MAO algorithm studied in [START_REF] Dimov | Monte Carlo Methods for Applied Scientists[END_REF][START_REF] Dimov | Monte Carlo Algorithms for Linear Problems[END_REF][START_REF] Dimov | Minimization of the Probable Error for Some Monte Carlo methods[END_REF][START_REF] Dimov | Robustness and applicability of Markov chain Monte Carlo algorithms for eigenvalue problems[END_REF]. Here we give a brief presentation of MAO. Suppose we have a Markov chain:

T = α 0 → α 1 → . . . α k → . . . , (7)
with n states. The random trajectory (chain) T k of length k starting in the state α 0 is defined as follows:

T k = α 0 → α 1 → . . . α j → . . . α k , (8)
where α j means the number of the state chosen, for j = 1, . . . , n.

Assume that

P (α 0 = α) = p α and P (α j = β|α j-1 = α) = p αβ , (9)
where p α is the probability that the chain starts in state α and p αβ is the transition probability to state β after being in state α. Probabilities p αβ define a transition matrix P . For the standard Markov chain Monte Carlo we normally require that

The above allows the construction of the following random variable, that is a unbiased estimator for the Problem P1:

θ[b] = b α 0 p 0 ∞ ∑ m=0 Q m b αm , (11)
where

Q 0 = 1; Q m = Q m-1 a α m-1 ,αm p α m-1 ,αm , m = 1, 2, . . . (12)
and α 0 , α 1 , . . . is a Markov chain on elements of the matrix A constructed by using an initial probability p 0 and a transition probability p α m-1 ,αm for choosing the element a α m-1 ,αm of the matrix A.

Now define the random variables Q m using the formula [START_REF] Dimov | Parallel computations of eigenvalues based on a Monte Carlo approach[END_REF]. One can see, that the random variables Q m , m = 1, . . . , i can also be considered as weights on the Markov chain [START_REF] Dimov | Monte Carlo Methods for Applied Scientists[END_REF].

One possible possible permissible densities can be chosen in the following way:

p = {p α } n α=1 ∈ P b , p α = |b α | ∥ b ∥ ; P = {p αβ } n α,β=1 ∈ P A , p αβ = |a αβ | ∥ a α ∥ , α = 1, . . . , n. (13)
Such a choice of the initial density vector and the transition density matrix leads to an Almost Optimal Monte Carlo (MAO) algorithm (see [START_REF] Dimov | Minimization of the Probable Error for Some Monte Carlo methods[END_REF][START_REF] Dimov | Monte Carlo Methods for Applied Scientists[END_REF]). The initial density vector p = {p α } n α=1 is called almost optimal initial density vector and the transition density matrix P = {p αβ } n α,β=1 is called almost optimal density matrix [START_REF] Dimov | Minimization of the Probable Error for Some Monte Carlo methods[END_REF].

Let us consider Monte Carlo algorithms with absorbing states: instead of the finite random trajectory T i in our algorithms we consider an infinite trajectory with a state coordinate δ i (i = 1, 2, . . .). Assume δ i = 0 if the trajectory is broken (absorbed) and δ i = 1 in other cases. Let

∆ i = δ 0 × δ 1 × . . . × δ i .
So, ∆ i = 1 up to the first break of the trajectory and ∆ q = 0 after that.

It is easy to show, that under the conditions (i) and (ii), the following equalities are fulfilled:

E{Q i f k i } = (h, A i b), i = 1, 2, . . . ; E{ N ∑ i=0 ∆ i Q i b k i } = (h, x), (P 1), E{ ∑ i|k i =r ′ ∆ i Q i } = c rr ′ , (P 3),
where (i|k i = r ′ ) means a summation only for weights Q i for which

k i = r ′ and C = {c rr ′ } n r,r ′ =1 .
The new WE algorithm studied in this paper was proposed by the second author. One important difference from the previously known Markov chain Monte Carlo algorithms is the special choice of the absorbtion probabilities.

In this case the Markov chain has n+1 states, {1, . . . , n, n+1}. The transition density matrix [START_REF] Dimov | Iterative Monte Carlo algorithms for linear algebra problems[END_REF] is such that ∑ n β=1 p αβ ≤ 1 and the probability for absorbtion at each stage accept the starting one is p α,n+1 = p α = 1 -

∑ n β=1 p αβ , α = 1, . . . , n.

Description of the probabilistic representation for the WE algorithm

We start the description of the probabilistic representation of the algorithm with two simple examples. The first example is for the case when all the coefficients of the matrix are non-negative. The second example deals with the case when some matrix coefficients may be negative real numbers. In this section we give a probabilistic representation for the WE algorithm for realvalued matrices, as well as for complex-valued ones.

A simple example: the positive case

We start with a simple example of positive case meaning that a ij ≥ 0, i, j = 1, . . . , n. For simplicity, we also assume that

∑ n j=1 a ij ≤ 1.
We describe on a very simple example how we obtain our probabilistic representation. We consider the system of two equations:

x 1 = 1 2 x 1 + 1 4 x 2 + 1, x 2 = 1 3 x 1 + 1 3 x 2 + 2
with unknowns x 1 and x 2 . We have obviously

x 1 = 1 + E(X)
where

P (X = x 1 ) = 1 2 , P (X = x 2 ) = 1 4 , P (X = 0) = 1 4
and

x 2 = 2 + E(Y )
where

P (Y = x 1 ) = 1 3 , P (Y = x 2 ) = 1 3 , P (Y = 0) = 1 3 .
We will use only one sample to approximate the expectations and compute a score initialized by zero along a random walk on the two equations. For instance, to approximate x 1 we add 1 to the score of our walk and then either we stop with probability 1 4 because we know the value of X which is zero or either we need again to approximate x(1) or u(2) with probability 1 2 or 1 4 . If the walk continues and that we need to approximate x(2), we add 2 to the score of our walk, we stop with probability 1 3 or either we need again to approximate x 1 or x(2) with probability 1 3 . The walk continues until one of the random variables X or Y takes the value zero and the score is incremented along the walk as previously. The score is an unbiased estimator of x 1 and we make the average of the scores of independent walks to obtain the Monte Carlo approximation of x 1 .

This algorithm is very similar to the ones used to compute the Feynman-Kac representations of boundary value problems by means of discrete or continuous walks like walk on discrete grids or the walk on spheres method. At each step of the walk, the solution at the current position x is equal to the expectation of the solution on a discrete or continuous set plus eventually a source term. The key ingredient is the double randomization principle which says that we can replace the expectation of the solution on this set by the solution at one random point picked according to the appropriate distribution. For example, for the walk on spheres methods this distribution is the uniform law on the sphere centered at x and of radius the distance to the boundary.

Simple example: general real matrix

We slightly modify our system by considering now the system of two equations

x 1 = 1 2 x 1 - 1 4 x 2 + 1, x 2 = 1 3 x 1 + 1 3 x 2 + 2
with unknowns x 1 and x 2 . We have x 1 = 1 + E(W ) where

P (W = x 1 ) = 1 2 , P (W = -x 2 ) = 1 4 , P (W = 0) = 1 4
and x 2 = 2 + E(Z) where

P (Z = x 1 ) = 1 3 , P (Z = x 2 ) = 1 3 , P (Z = 0) = 1 3 .
To approximate x 1 we add 1 to the score of our walk and then either we stop with probability 1 4 because we know the value of X which is zero or either we need to approximate x 1 or -x 2 with probability 1 2 or 1 4 . If the walk continues and that we need to approximate -x 2 , we add -2 to the score of our walk, we stop with probability 1 3 or either we need to approximate -x 1 or -x 2 with probability 1 3 . If the walk continues again and that we need to approximate -x 1 , we add -1 to the score of our walk, we stop with probability 1 4 or either we need to approximate -x 1 or x 2 with probability 1 2 or 1 4 . The walk continues until one of the random variables W or Z takes the value zero and the score is incremented along the walk as previously. The same kind of reasoning holds for a complex-valued matrix in the representation given later on in Subsection 3.4.

Probabilistic representation for real matrices

Consider a real linear system of the form u = Au + b where the matrix A of size n is such that ϱ(A) < 1, its coefficients a i,j are real numbers and

n ∑ j=1 |a i,j | ≤ 1, ∀1 ≤ i ≤ n.
We now define a Markov chain T k with n + 1 states α 1 , . . . , α n , n + 1, such that

P (α k+1 = j/α k = i) = |a i,j | if i ̸ = n + 1 and P (α k+1 = n + 1/α k = n + 1) = 1.
We also define a vector c such that c(i) = b(i) if 1 ≤ i ≤ n and c(n + 1) = 0. Denote by τ = (α 0 , α 1 , . . . , α k , n + 1) a random trajectory that starts at the initial state α 0 < n + 1 and passes through (α 1 , . . . , α k ) until the absorbing state α k+1 = n + 1. The probability to follow the trajectory τ is P (τ ) = p α 0 p α 0 α 1 , . . . p α k-1,k α k p α k . We use the MAO algorithm defined by ( 13) for the initial density vector p = {p α } n α=1 and for the transition density matrix P = {p αβ } n α,β=1 , as well. The weights Q α are defined in the same way as in [START_REF] Dimov | Parallel computations of eigenvalues based on a Monte Carlo approach[END_REF], i.e.,

Q m = Q m-1 a α m-1 ,αm p α m-1 ,αm , m = 1, . . . , k, Q 0 = c α 0 p α 0 . ( 14
)
The estimator θ α (τ ) can be presented as

θ α (τ ) = c α + Q k aα k α pα k , α = 1, . . . n taken with a probability P (τ ) = p α 0 p α 0 α 1 , . . . p α k-1,k α k p α k .
Now, we can prove the following theorem.

Theorem 3.1 The random variable θ α (τ ) is an unbiased estimator of x α , i.e. E{θ α (τ )} = x α . P r o o f. E{θ α (τ )}= ∑ τ =(α 0 ,...,α k ,n+1) θ α (τ )P (τ ) = ∑ τ =(α 0 ,...,α k ,n+1) ( c α + Q k (τ ) a α k α p α k ) P (τ ) = c α + ∑ τ =(α 0 ,...,α k ,n+1) Q k (τ ) a α k α p α k P (τ ) = c α + ∑ τ =(α 0 ,...,α k ,n+1) c α 0 p α 0 a α 0 α 1 . . . a α k-1 α k p α 0 α 1 . . . p α k-1 α k a α k α p α k p α 0 × p α 0 α 1 . . . p α k-1 α k p α k = c α + ∞ ∑ k=0 n ∑ α 0 =1 . . . n ∑ α k =1 a α k α a α k-1 α k . . . a α 0 α 1 c α 0 = c α + (Ac) α + (A 2 c) α + . . . = ( ∞ ∑ i=0 ( A i c ) ) α . (15)
Since ϱ(A) < 1, the last sum in ( 15) is finite and

( ∞ ∑ i=0 ( A i c ) ) α = x α , α = 1, . . . , n.
Thus, E{θ α (τ )} = x α . ♢

Let as stress on the fact that among elements c α 0 , a α 0 α 1 . . . a α k-1 α k there are no elements equal to zero because of the special choice of permissible distributions p α and p αβ defined by MAO densities [START_REF] Dimov | A Power Method with Monte Carlo Iterations[END_REF]. The rules [START_REF] Dimov | A Power Method with Monte Carlo Iterations[END_REF] ensure that the Markov chain visits non-zero elements only.

One may also observe that the conjugate to

θ α (τ ) random variable θ * α (τ ) = ∑ k i=0 Q i aα i α pα i c α
P r o o f. Equality ( 16) can easily be proved using properties of the mean value.

Using the variance properties one can also show that ( 17) is true. ♢ Moreover, using the Kolmogorov theorem and the fact that the sequences (θ α,N ) N ≥1 are independent and also independently distributed with a finite mathematical expectation equal to x α one can show that θα,N converges almost surely to x α :

P ( lim N →∞ θα,N = x α ) = 1.
The theorems proved (Th.3.1 and Th.3.2) with the last fact allows us to use the Monte Carlo method with the random variable θ α for solving algebraic systems and ensure that the estimate is unbias and converges almost surely to the solution. These facts do not give us a measure of the rate of convergence and do not answer the question how one could control the quality of the Monte Carlo algorithm? At the same time it is clear that the quality of the algorithm, meaning how fast is the convergence and how small is the probability error, should depend on the properties of the matrix A, and possibly of the righthand side vector b. Our study shows that the balancing of the matrix together with the matrix spectral radius is an important issue for the quality of the Monte Carlo algorithm. And this is important for all three problems under consideration. To analyse this issue let us consider the variance of the random variable θ α (τ ) for solving Problem P1, i.e., evaluating the inner product

J(x) = (v, x) = n ∑ i=1 v i x i of a given vector v = (v 1 , . . . , v n ) t ∈ IR n×1
and the solution x of the system (2):

x = Ax + b, A ∈ IR n×n , b = (b 1 , . . . , b n ) t ∈ IR n×1
. By chousing v to be an unit vector e (j) ≡ (0, . . . , 0, 1 j , 0, . . . , 0) all elements, of which are zeros except the j th element e (j) j , which is equal to 1 the functional J(x) becomes equal to x j . Denote by Ā the matrix containing the absolute values of elements of a given matrix A: Ā = {|a ij |} n i,j=1 and by ĉ the vector of squared entrances of the vector c, i.e., ĉ = {c 2 i } n+1 i=1 . The pair of MAO density distributions (13) define a finite chain of vector and matrix entrances:

c α 0 → a α 0 α 1 → . . . → a α k-1 α k . ( 18
)
The latter chain induces the following product of matrix/vector entrances and norms:

A k c = c α 0 k ∏ s=1 a α s-1 αs ,
where the vector c ∈ IR n×1 is defined as

c(i) = b(i) if 1 ≤ i ≤ n and c(n + 1) = 0.
Obviously, the variance of the random variable ψ k α (τ ) defined as

ψ k α (τ ) = c α 0 p α 0 a α 0 α 1 p α 0 α 1 a α 1 α 2 p α 1 α 2 . . . a α k-1 α k p α k-1 α k c α k p α k = A k c c α k P k (τ )
plays an important role in defining the quality of the proposed WE Monte Carlo algorithm. Smaller is the variance of D{ψ k α (τ )}, better convergency of the algorithm will be achieved.

The next theorem gives the structure of the variance for the WE algorithm.

Theorem 3.3 D{ψ k α (τ )} = c α 0 p α 0 p α ( Āk c ĉ) α -(A k c c) 2 α . P r o o f.
We deal with the variance

D{ψ k α (τ )} = E{(ψ k α (τ )) 2 } -(E{ψ k α (τ )) 2 . ( 19
)
Consider the first term of [START_REF] Halton | Sequential Monte Carlo[END_REF].

E{(ψ k α (τ )) 2 } = E { c 2 α 0 p 2 α 0 a 2 α 0 α 1 . . . a 2 α k-1 α k p 2 α 0 α 1 . . . p 2 α k-1 α k c 2 α k p α k } = n ∑ α 0 ,...α k =1 c 2 α 0 p 2 α 0 a 2 α 0 α 1 . . . a 2 α k-1 α k p 2 α 0 α 1 . . . p 2 α k-1 α k c 2 α k p 2 α k p α 0 p α 0 α 1 . . . p α k-1 α k p α k = n ∑ α 0 ,...α k =1 c 2 α 0 p α 0 |a α 0 α 1 | . . . |a α k-1 α k | c 2 α k p α k = n ∑ α 0 ,...α k =1 |c α 0 | p α 0 p α k |c α 0 ||a α 0 α 1 | . . . |a α k-1 α k |c 2 α k = |c α 0 | p α 0 p α × ( Āk c ĉ) α .
In the same way one can show that the second term of ( 19) is equal to (A k c c) 2 α .

Thus,

D{ψ k α (τ )} = cα 0 pα 0 pα ( Āk c ĉ) α -(A k c c) 2 α ♢
The problem of robustness of Monte Carlo algorithms is considered in [START_REF] Dimov | Robustness and applicability of Markov chain Monte Carlo algorithms for eigenvalue problems[END_REF] and [START_REF] Dimov | Monte Carlo Methods for Applied Scientists[END_REF]. We will not go into details here, but we should mention that according to the definition given in [START_REF] Dimov | Monte Carlo Methods for Applied Scientists[END_REF] the proposed WE algorithm is robust. In [START_REF] Dimov | Monte Carlo Methods for Applied Scientists[END_REF] the problem of existing of interpolation MC algorithms is also considered. Simply specking, MC algorithm for which the probability error is zero is called interpolation MC algorithm. Now we can formulate an important corollary that gives a sufficient condition for constructing an interpolation MC algorithm.

Corollary 3.1 Consider a perfectly balanced stochastic matrix

A = {a ij } n i,j=1
, where a ij = 1/n, for i, j = 1, . . . , n and the vector c = (1, . . . , 1) t . Then MC algorithm defined by density distribution P k (τ ) is an interpolation MC algorithm.

P r o o f. To prove the corollary it is sufficient to show that the variance D{ψ k α (τ )} is zero. Obviously,

Ac =        1 n . . . 1 n . . . 1 n . . . 1 n               1 . . . 1        =        1 . . . 1       
, and also,

A k c =        1 . . . 1        . Since Ā =        1 n . . . 1 n . . . 1 n . . . 1 n        = A and ĉ = {c 2 i } n i=1 = c, we have: |cα 0 | pα 0 pα × ( Āk c ĉ) α = 1. Thus, we proved that D{ψ k α (τ )} = 0 ♢
First of all, let us mention that solving systems with perfectly balanced stochastic matrices does not make sense. Nevertheless, it is possible to consider matrix-vector iterations A k c since they are the basics for the von Neumann series that approximate the solution of systems of linear algebraic equations under some conditions given in Section 2. The idea of this consideration was to demonstrate that, as closer the matrix A is to the perfectly balanced matrix, as smaller is the probability error of the algorithm.

Probabilistic representation for complex-valued matrices

Assume that one wants to solve a linear system of the form u = Au + b, where the square matrix A ∈ C n×n with complex-valued entrances is of size n. Assume also that ϱ(A) < 1, its coefficients a i,j verify

∑ n j=1 |a i,j | ≤ 1, ∀1 ≤ i ≤ n.
We also define a Markov T k with n + 1 states α 1 , . . . , α n , n + 1, such that

P (α k+1 = j/α k = i) = |a i,j |
where f j ≡ e j ≡ (0, . . . , 0, 1 j , 0, . . . , 0) and g j ≡ (g j1 , g j2 , . . . , g jn ) t is the j-th column of the inverse matrix G = B -1 .

Here we deal with the matrix A = {a ij } n ij=1 , such that

A = I -DB, ( 22
)
where D is a diagonal matrix D = diag(d 1 , . . . , d n ) and d i = γ b ii , i = 1, . . . , n, and γ ∈ (0, 1] is a parameter that can be used to accelerate the convergence. The system [START_REF] Halton | Sequential Monte Carlo[END_REF] can be presented in the following form:

u = Au + b, (23)
where B = Df.

Let us suppose that the matrix B has diagonally dominant property. In fact, this condition is too strong and the presented algorithms work for more general matrices, as it will be shown later on. Obviously, if B is a diagonally dominant matrix, then the elements of the matrix A must satisfy the following condition:

n ∑ j=1 |a ij | ≤ 1 i = 1, . . . , n. (24)

A new WE Monte Carlo algorithm for computing one component of the solution

Here we describe a new WE Monte Carlo algorithm for computing the component u i 0 of the solution.

Algorithm 4.1 :

Computing one component x i 0 of the solution x i , i = 1, . . . n 1. Initialization Input initial data: the matrix B, the vector f , the constant γ and the number of random trajectories N .

Preliminary calculations (preprocessing):

2.1. Compute the matrix A using the parameter γ ∈ (0, 1]:

{a ij } n i,j=1 =      1 -γ when i = j -γ b ij b ii when i ̸ = j .

2.2.

Compute the vector asum:

asum(i) = n ∑ j=1 |a ij | for i = 1, 2, . . . , n.

Compute the transition probability matrix

P = {p ij } n i,j=1
, where

p ij = |a ij | asum(i)
, i = 1, 2, . . . , n j = 1, 2, . . . , n .

A new WE Monte Carlo algorithm for computing all components the solution

Here we describe the WE algorithm for computing all the components of the solution. For simplicity we will consider the case when all the matrix entrances are nonnegative. So, we assume that a ij ≥ 0, for i, j = 1, . . . n. We will skip the part of initialization and preprocessing, which is the same as in the previous subsection.

Algorithm 4.2 :

Computing all components x i , i = 1, . . . n of the solution 1. Initialization.

Preprocessing. 3. for

i = 1 to n do 3.1. S(i) := 0; V (i) := 0 3.2. for k = 1 to N do (MC loop) 3.1 set m := rand(1 : n) .
3.2 set test := 0; m 1 := 0; 3.3 while (test ̸ = 1 do:

3.3.1. V (m) := V (m) + 1; m = m 1 + 1; l(m 1 ) = m; 3.3.2.
for q = 1, m 1 do: S(l(q)) := S(l(q)) + b(l(q)); 3.3.3. if r > asum(m), then test = 1;

else chose j according to the probability density matrix p i,j ; 3.3.4. set m = j. 3.3.5. endif 3.4. endwhile 3.5. enddo 3.6. for j = 1, n do:

V (j) := max{1, V (j)}; S(j) = S(j) V (j) .

Sequential Monte Carlo

The sequential Monte Carlo method for linear systems has been introduced by John Halton in the sixties [START_REF] Halton | Sequential Monte Carlo[END_REF], and further developed in more recent works [START_REF] Halton | Sequential Monte Carlo[END_REF][START_REF] Halton | Monte Carlo integration with sequential stratification[END_REF][START_REF] Halton | Sequential Monte Carlo for linear systems -a practical summary[END_REF]. The main idea of the SMC is very simple. Assume that one needs to solve the problem:

x = Ax + b, A ∈ IR n×n , b = (b 1 , . . . , b n ) t ∈ IR n×1 .
We put: x = y +z, where y is a current estimate of the solution x using a given Monte Carlo algorithm, say WE algorithm and z is the correction. Then the correction, z satisfies z = Az +d, where d = Ay -y +b. At each stage, one uses a given Monte Carlo to estimate z, and so, the new estimate y. If the sequential computation of d is itself approximated, numerically or stochastically, then the expected time for this process to reach a given accuracy is linear per number of steps, but the number of steps is dramatically reduced.

Normally, the SMC leads to dramatic reduction of the computational time and to high accuracy. Some numerical results demonstrating this fact will be shown in the next section. The algorithms of this kind can successfully compete with optimal deterministic algorithms based on Krilov subspace methods and the conjugate gradient methods (CGM) (many modifications of these methods are available). But one should be careful, if the matrix-vector multiplication Ay is performed in a deterministic way, it needs O(n 2 ) operations. It means that for a high quality WE algorithm combined with SMC one needs to have very few sequential steps. In the next section some numerical results with high accuracy obtained in just several sequential steps will be presented.

Numerical tests and discussion

In this section we show numerical tests divided into several groups. In the first group we compare the results obtained by the proposed WE algorithm combined with SMC to the results obtained by standard deterministic Jacobi algorithm for some test matrices. In the second group of numerical experiments we study how the quality of the WE algorithm depends on the balancing of the matrices. In the first two groups of experiments we deal with matrices of different size n ranging from n = 10 to n = 25000. The matrices are either randomly generated, or created with some prescribed properties we want to have, in order to study the dependance of the convergency of the algorithm to different properties. In the third group of experiments we run our algorithm to some large arbitrary matrices describing real-life problems taken from some available collections of matrices. In this particular study we use matrices from the Harwell-Boeing Collection. It is a well-know fact that Monte Carlo methods are very efficient to deal with large and very large matrices. This is why we pay special attention to these matrices.

Comparison of WE with deterministic algorithms

In this subsection we compare the numerical results for solving systems of linear equations by WE Monte Carlo combined with SMC with the well known Jacobi method. By both methods we compute all the component of the solution. Here we deal with dense matrices of both types: randomly generated, as well as deterministic matrices. The randomly generated matrices are such that max i

∑ n j=1 |a ij | = max ∥ a i ∥= a ≤ 1.
In most of the graphs presented here the parameter a was chosen to be equal to 0.9. On the graphs we present the relative error of the solution of the first component x 1 , since the results for all other components are very similar. So, on the graphs we present the value of

ε(x i ) = x i,comp -x i,exact x i,exact , i = 1, . . . , n,
where x i,comp is the computed value of the i-th component of the solution vector, and x i,exact is the i-th component of the exact solution-vector. On Figure 1 results for the relative error of the solution for the proposed algorithm and for the Jacobi method are presented.

One can observe that for the proposed WE Monte Carlo method the convergence is much faster and just after 30 iterations the solution is very accurate: the relative error is about 10 -12 . At the same time, the convergence of Jacobi method is slower, and after 100 iterations the relative error is between 10 -2 and 10 -3 . Similar results are obtained for much large dense matrix of size In order to see if the residual can be considered as a good measure for the error we compare the relative error with the residual for a matrix of size n = 5000 (see. Figure 3). One can observe that the curves presenting the relative error and the residual are very close to each other.

Role of balancing for the quality of the algorithms

Let us mention that the balancing of matrices is a vary important issue for the deterministic algorithms. That is why for any arbitrary matrix a balancing procedure is needed as a preprocessing. Monte Carlo algorithms are not that sensitive to the balancing, but nevertheless, as it was shown in Subsection 3.3 the balancing affects the convergence of the WE Monte Carlo algorithm. The first example is shown on Figure 4 where we consider two matrices of small size n = 10. The balanced matrix is such that a ij = 0.1, such that the parameter a defined above is equal to 0.9. For the unbalanced matrix the parameter a is equal to 0.88, but the unbalanced matrix is extremity unbalanced (the ration between the largest and the smallest by modulo elements is equal to 40). Let us stress on the fact that this ratio is considered large when it is equal to 5.

The unbalanced matrix we deal with is the following: Obviously, for the unbalanced matrix there is no convergence at all. For the balanced matrix the convergence is fast, but the red curve is rough. The reason for that is that we use a very small number of random trajectories (N = 5n). Now we can consider the same matrices, but let us increase the number of the random trajectories to N = 500. At the same time we may decrease the number of the random trajectories (and also the computational time) from N = 50 to N = 30 for the balanced matrix. The results of this numerical experiments are presented on Figure 5.

One can clearly see, that the WE algorithm starts to converge and after 35 to 40 iterations we have 3 correct digits. For the balanced matrix the convergence is still very fast, but again the red curve is rough, because of the small number of random trajectories. For a randomly generated matrices of size n = 5000 similar results are presented on Figure 6. The unbalanced matrix is generated by perturbing randomly the elements of the balanced matrix. The green and the blue curves correspond to the same unbalanced matrix, but the numerical results presented by the blue curve are obtained by using 10 times more random trajectories. In such a way, one can see that increasing the number N of random trajectories one can improve the convergence, as it is clear from the theory. In our particular numerical experiments shown on Figure 6 the green curve is obtained after a number of trajectories N = 5n, while the blue curve is based on N = 50n number of trajectories. One can see that in the first case the algorithm does not converge, while in the second case it converges. These experiments show that the WE Monte Carlo has another degree of freedom (which does not exist in deterministic algorithms) to improve the convergence even for very bad balanced matrices. The price we pay to improve the convergence is the increased complexity (or, the computational time) since we increase the number of simulation trajectories. It is interesting to see how the disbanding affects the convergence of the WE Monte Carlo. We have performed the following numerical experiment: taking randomly generated matrices of differen size we were perturbing the matrix elements randomly by some percentage making the matrix more and more imbalanced. On Figure 7 we present the results for the matrix of size n = 100. The randomly generated entrances of the matrix are additionally perturbed by 5% and by 20%.

Comparison of WE with the preconditioned conjugant gradient (PCG) method

We compare our results with the optimal deterministic preconditioned conjugant gradient (PCG) method [START_REF] Golub | Inexact Preconditioned Conjugate Gradient Method with Inner-Outer Iteration[END_REF][START_REF] Knyazev | Steepest Descent and Conjugate Gradient Methods with Variable Preconditioning[END_REF][START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. Here is more precipe information about our implementation of the PCG method. We want to solve the linear system of equations Bx = b by means of the PCG iterative method. The input arguments are the matrix B, which in our case is the square large matrix N OS4 taken from the well-known Harwell-Boeing Collection. B should be symmetric and positive definite. In our implementation we use the following parameters:

• b is the right-hand side vector.

• tol is the required relative tolerance for the residual error, namely b -Bx.

The iteration stops if ∥ b -Bx ∥≤ tol * ∥ b ∥. • maxit is the maximum allowable number of iterations.

• m = m1 * m2 is the (left) preconditioning matrix, so that the iteration is (theoretically) equivalent to solving by PCG P x = m b, with P = m B. • x 0 is the initial guess.

A comparison for convergency of the WE Monte Carlo and the PCG for the matrix N OS4 is presented in the Figure 8. This particular matrix is taken from an application connected to finite element approximation of a problem describing a beam structure in constructive mechanics [START_REF]Matrix market, NOS4: Lanczos with partial reorthogonalization[END_REF]. The results presented on Figure 8 show that the convergence for the WE Monte Carlo is better: the curve presenting the residual with the number of iterations is smoother and goes down to 10 -6 -10 -7 , while the curve presenting the PCG achieves an accuracy of 10 -3 . Let us stress on the fact that for the experiments presented the needed accuracy is set to tol = 10 -8 . One can not guarantee that such an effect happens for every matrices, but there are cases in which the WE Monte Carlo performs better that the PCG.

Concluding remarks

A new Walk on Equations Monte Carlo algorithm for solving linear algebra problems is presented and studied. The algorithm can be used for evaluating all the components of the solution of systems. The algorithm can also be used to approximate the matrix inversion of square matrices. The WE Monte Carlo is applied to real-valued, as well as to complex-valued matrices. It is proved that:

• the proposed random variable based on random walk on equations is an unbiased estimate of the solution vector; • for the Monte Carlo estimator, the mean value of N realizations of the random variable, it is proven that the variance vanishes when N goes to infinity and converges almost surely to the solution; • the structure of the relative stochastic error is analysed; it is shown that an interpolation Monte Carlo algorithm is possible for perfectly balanced matrices.

The analysis of the structure of the stochastic error shows that the rate of convergence of the proposed WE Monte Carlo algorithm depend on the balancing of the matrix. This dependance is not that strong as for deterministic algorithms, but still exists.

A number of numerical experiments are performed. The analysis of the results show that:

• the proposed WE Monte Carlo algorithm combined with the sequential Monte Carlo for computing all the components of the solution converges much faster than some well-known deterministic iterative methods, like Jacobi method; this is true for matrices of different size and the effect is much bigger for larger matrices of size up to n = 25000; • the balancing of matrices is still important for the convergency; many numerical results demonstrate that one should pay special attention to the balancing; • the comparison of numerical results obtained by the WE Monte Carlo and the best optimal deterministic method, the preconditioned conjugant gradient, show that for some matrices WE Monte Carlo gives better results.

p

  αβ = 1 for any α = 1, . . . , n.

  has the same mean value, i.e., x α . It can be immediately checked following the same technique as in the proof of Theorem 3.1. Assume one can compute N values of the random variable θ α , namely θ α,i , i = 1, . . . , N . Let us consider the value θα,N = 1 N ∑ N i=1 θ α,i as a Monte Carlo approximation of the component x α of the solution. For the random variable θα,N the following theorem is valid: Theorem 3.2 The random variable θα,N is an unbiased estimator of x α , i.e. E{ θα,N } = x α (16) and the variance D{ θα,N } vanishes when N goes to infinity, i.e lim N →∞ D{ θα,N } = 0.

3 .

 3 Set S := 0 . 4. for k = 1 to N do (MC loop) 4.1 set m := i 0 . 4.2 set S := S + b(m). 4.3. test = 0; sign = 1 4.4 while (test ̸ = 0 do: 4.4.1. generate an uniformly distributed random variable r ∈ (0, 1) 4.4.2. if r ≥ asum(m) then test = 1; 4.4.3. else chose the index j according to the probability density matrix p ij ; 4.4.4. set sign = sign * sign{a mj }; m = j; 4.4.5 update S := S + sign * b m ; 4.4.6. endif . 4.5 endwhile 5. enddo

Fig. 1 .Fig. 2 .

 12 Fig.1. Relative error for the WE Monte Carlo for all components of the solution and Jacobi method for a matrix of size 5000 × 5000. The number of random trajectories for MC is 5n.

Fig. 3 .

 3 Fig. 3. Comparison the relative error with the residual for the WE Monte Carlo 100 × 100. The number of random trajectories for MC is 5n.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Residual for the WE Monte Carlo for balanced and extremity unbalanced matrix of size n = 10. The number of random trajectories is N = 5n.

Fig. 6 .

 6 Fig.6. Residual for the WE Monte Carlo; the matrix of size is 5000 × 5000. The number of random trajectories for MC is 5n.

Fig. 7 .

 7 Fig. 7. Residual for the WE Monte Carlo; the matrix of size is 100 × 100. Matrix elements are additionally perturbed by 5% and by 20%.

Fig. 8 .

 8 Fig. 8. Comparison of the WE Monte Carlo and the PCG method for a matrix N OS4 from the Harwell-Boeing Collection.

  Consider also a transition density matrix P = {p ij } n i,j=1 ∈ IR n×n , such that p ij ≥ 0, i, j = 1, . . . , n and Define sets of permissible densities P b and P A .

	∑ n j=1 p ij = 1, for
	any i = 1, . . . , n. Definition 2.1 The initial density vector p = {p i } n i=1 is called permissible to the vector b = {b i } n i=1 ∈ IR n , i.e. p ∈ P b , if
	
	 
	 

Acknowledgment

The research reported in this paper is partly supported by the European Commission under FP7 project AComIn (FP7 REGPOT-2012-2013-1, by the Bulgarian NSF Grant DCVP 02/1, as well as by the Université du Sud Toulon-Var.

⋆ This work was partially supported by the European Commission under FP7

project AComIn (FP7 REGPOT-2012-2013-1, by the Bulgarian NSF Grant DCVP 02/1, as well as by the Université du Sud Toulon-Var.

if i ̸ = n + 1 and

Like in the real case, we also define a vector c such that c(i) = b(i) if 1 ≤ i ≤ n and c(n + 1) = 0. Denote by τ = (α 0 , α 1 , . . . , α k , n + 1) a random trajectory that starts at the initial state α 0 < n + 1 and passes through (α 1 , . . . , α k ) until the absorbing state α k+1 = n + 1. The probability to follow the trajectory τ is

Additionally to the Markov chain T k and to the vector c we define a random variable S α k such that

in the real case. Then, we have the following probabilistic representation:

The proof of this probabilistic representation is similar to the proof in the real case given in Theorem 3.1.

Description of the walk on equations (WE) algorithms

We start this section with a remark of how the matrix inversion problem (see Problem P3 presented in Section 2) can be treated by the similar algorithm for linear systems of equations. Consider the following system of linear equations:

The inverse matrix problem is equivalent to solving n-times the problem (20), i.e.

Bg j = f j , j = 1, . . . , n


Jan Christopher Bernhammer

Gilles Frison

Han Vinh Huynh

Pincer versus Pseudo-Pincer: Isomerism in Palladium(II) Complexes bearing κ 3 C,S,C ligands

In NHC pincer complexes incorporatig a hemilabile donor site, there exists an equilibrium between the true pincer form and a pseudopincer coordination isomer. The influence of the NHC moieties on this isomerism has been studied by DFT calculations.

Tridentate, pincer and pincer-type ligands have gained increasing popularity due to the robustness and thermal stability the rigid, κ 3 -mer or κ 3 -(pseudo)-mer coordination geometry these ligands confer to their complexes. Transition metal complexes incorporation pincer ligands have been shown to be versatile and useful catalysts in organic synthesis. [START_REF]The Chemistry of Pincer Compounds[END_REF] Beyond catalysis, pincer complexes have also been utilized successfully for the development of new materials. 2 The fine-tuning of steric and electronic properties of pincer ligands -be it by modifications in the backbone or changes to 20 the nature of the donor sites -has attracted considerable research interest, in order to provide coordination spheres for each and every application. 3 While catalysis often involves the coordinations sites on the metal centre not occupied by the pincer ligand, there are also cases in which partial decoordination of a hemilabile donor takes place. The advantage is the rapid creation of a coordinatively and electronically unsaturated metal centre, open for incoming substrates, while retaining the tethered donor functionality in the vicinity in order to allow for rapid recoordination and stabilisation of the active species. As a result, highly active, yet stable catalysts are obtained. 4 Frequently, the donor functionalities in the sidechains of the pincer ligands are the hemilabile sites, 5 but cases in which the central donor moiety exhibits hemilabile behaviour have been reported as well. 6 However, to control hemilabile behaviour, the governing geometric and electronic factors still need to be understood. 7 N-heterocyclic carbenes (NHCs) have become ubiquitously employed ligands for transition metal catalysts due to their strong σ-basicity, the robustness of their complexes, and the 40 ease with which they can be synthesized. 8 Their steric bulk and electronic properties can be fine-tuned rapidly by modifying the backbone and side chains, 9 and additional donor functionalities can be introduced in the side chains, which gives rise to chelating and pincer-type ligands. Donor 45 moieties based on heteroatoms such as nitrogen, oxygen, phosphorous, and sulfur have been reported. 10,11 Especially the soft thioether functionalities are know to exhibit hemilabile behavior in their complexes. 12 It has been shown that in the presence of coordinating 50 counteranions, complexes of CSC pincer ligands are isolated either as κ 2 -trans-C,C pseudo-pincer complexes or in their κ 3mer-C,S,C pincer form. 13 The electron-donating abilities of the carbene moieties were suggested to have a pivotal influence on this isomerism. To study this electronic effect, 55 we examined the electronic properties of NHCs with markedly different donor strenghts and the preference for either the pincer or the pseudo-pincer form by means of computational methods.

The relative Gibbs free energies ΔG of a neutral 60 pseudopincer species A and the corresponding cationic pincer C and a bromide ion are difficult to obtain due to the significant impact of solvatation on the ion pair. Instead, the Gibbs free energy ΔG R of the reaction of a pincer and pseudopincer complex was used as an indication of preference 65 for the pincer form (Scheme 1). During this fictitious reaction, the number of ionic species stays constant throughout the reaction, and solvatation effects as well as the Gibbs free energy of the bromide counteranion can be neglected as they are almost identical for starting materials and products.

70

The resulting scale of pincer preference energy is arbitrary in the sense that the absolute values will depend on the choice of reference system. We decided to use complexes B and D featuring dichloroimidazolin-2-ylidene moieties as probe. This system was the most weakly donating ligand under scrutiny, and the respective complexes are known to prefer the pseudopincer form. A negative value for ΔG R will thus signify a higher preference for the pincer form in comparison to the reference system. Both the pincer and pseudopincer forms of all examined complexes were optimised at the B3LYP/cc-pVDZ level of theory, 14,15,16 and their respective Gibbs free energies calculated at the same level of theory (see Supplementary Information for computational details).

To avoid the influence of variations in steric bulk 10 associated with changing the substituents at nitrogen and to limit computational cost, NHC side chains were simplified to sterically unassuming methyl groups, and the effect of five vastly different NHC backbones (cf. figure 2) was examined, two of which have already been studied experimentally. The optimised structures are comparable to the reported molecular structures. 13 All complexes have a distorted square 20 planar coordination geometry of the palladium center (Figure 1). In the pincer complexes, the metallacycles adopt a distorted boat conformation and the planes defined by the NHC rings are almost coplanar, while they showed a twist of 26-43° in the corresponding pseudopincer forms. From these optimized structures, ΔG R can be obtained.

No experimental donor strength values are available for the ligands in the calculated complexes. However, it has been shown that the Kohn-Sham eigenvalue of the σ-lone pair orbital ε(σ-HOMO) correlates well with the donor strength. 17 To limit computational cost, orbital energies were calculated for truncated ligands at the B3LYP/auc-cc-pVTZ//B3LYP/cc-pVDZ level (Figure 2). 14,16,18 In all cases, the highest orbital with σ-symmetry was found to be the lone pair at C carbene (Figure 3). Orbital energies for the σ-lone pair orbital were found to 40 range from -6.56 to -5.52 eV, and confirmed the expected increase in donor strength from 1 to 5 (Table 1). The preference for the pincer form, represented by ΔG R of the test reaction, increased in the same order, and Figure 4 shows the strong correlation between donor strength and ΔG R (R 2 = 45 0.9931). It is apparent that more strongly donating carbene moieties in the CSC ligand favour the pincer form. The cationic pincer 50 complexes are better stabilised by stronger NHCs than the neutral pseudopincer form, and this difference in stabilization increases with the donor strength. Additionally, entropy and solvation also favour the formation of the ionic pincer over the pseudopincer form. On the other hand, the neutral 55 pseudopincer complexes are preferred if the NHC moieties transfer less electron density towards the metal atom due to the fact that bromide is a better ligand than the thioether. In a more formalistic approach, the [PdBr(NHC) 2 ] + fragment can be understood as a Lewis acid interacting with a either the thioether moiety or a bromido ligand as a Lewis base. The internal charge redistribution in the Lewis acid upon charge transfer from the Lewis base is energetically more favourable with more weakly donating NHC moieties, which allows the coordination of the more strongly donating bromido ligand. 19 -6,6 -6,4 -6,2 -6,0 -5,8 -5,6 -5,4 - In summary, the relationship between donor strength and pincer/pseudo-pincer isomerism has been unambigously demonstrated in good agreement with experimental observations. Calculated or experimentally determined ε(σ-HOMO) values can thus serve as indication for the hemilability of the central donor site in similar complexes, allowing for a more rational ligand design and the synthesis of taylor-made pincer complexes.

We thank the vice-presidency for external relations in École polytechnique for support from the international 20 internship program and the National University of Singapore and the Singapore Ministry of Education for financial support (WBS R-143-000-410-112) and a SINGA scholarship for J.C.B. This work was performed using HPC resources GENCI-CINES (Grant 2013-x086894).

25

Notes and references

a Department of Chemistry, National University of Singapore, 3 Science Drive 3, 117543, Singapore. Fax: +65 6779 1691; Tel: +65 6516

Scheme 1

 1 Scheme 1 Homodesmotic pincer/pseudo-pincer exchange reaction

13

 13

Fig. 1

 1 Fig. 1 Optimised pseudo-pincer and pincer geometries based on the dichloromidazolin-2-ylidene moiety

Fig. 2

 2 Fig. 2 Truncated ligands with different NHC moieties

  This journal is © The Royal Society of Chemistry [year] Journal Name, [year], [vol], 00-00 | 3

Fig. 4

 4 Fig. 4 Correlation ε(σ-HOMO) [eV] and ΔGR [kJ/mol]

Table 1 .

 1 ε(σ-HOMO) of the truncated ligand and ΔGR

	NHC moiety	ε(σ-HOMO) [eV]	ΔGR [kJ/mol]
	Cl2-imidazole	-6.56	0.0
	1,2,4-triazole	-6.52	-5.4
	Imidazole	-6.12	-26.5
	Imidazoline	-5.94	-32.8
	Pyrazole	-5.52	-57.0

  2670; E-mail:chmhhv @nus.edu.sg b Laboratoire de Chimie Moleculaire, Department of Chemistry, Ecole

	polytechnique and CNRS, 91128 Palaiseau Cedex, France. Fax: +33 1
	6933 4803; Tel: +33 1 6933 4834; E-mail:
	gilles.frison@polytechnique.org
	† Electronic Supplementary Information (ESI) available. See
	DOI: 10.1039/b000000x/