NC solving of a system of linear ordinary differential equations in several unknowns

D. Grigoriev

An NC algorithm is described for reducing a system of linear ordinary differential equations in several unknowns to the standard basis form

Introduction

We consider a problem of solvability of a system of linear ordinary differential equations in several unknowns

∑_{j} L_{ij}u_{j} = b_{i}

where b_{i} ∈ ℚ(X) and L_{ij} = ∑_{k} f_{k} \frac{d^{k}}{dx^{k}} are linear ordinary differential operators with the rational coefficients f_{k} ∈ ℚ(X). We consider a solvability of the system in the unknowns u_{j} in the differential closure of ℂ(X) (in fact, as we deal with the linear operators it is equivalent to the solvability in Picard-Vessiot closure of ℂ(X)) (see [K] ), in which any (resp. linear) ordinary differential equation has a solution. In other words, solvability in Picard-Vessiot closure means that the system cannot be brought to a contradiction by equivalent transformations over the ring R = ℂ(X) [\frac{d}{dX} ] of the linear differential operators, or more precisely, that the ideal in the ring of differential polynomials in {u_{j}}, generated by the differential polynomials {∑_{j} L_{ij}u_{j} − b_{i}}, differs from the unit one.

Remark that this problem is a particular case of the problem of solvability (over differential closure) of a system of non-linear ordinary differential equations in several unknowns (more general, a quantifier elimination problem for these systems), for which an algorithm with elementary complexity (more precisely, double-exponential) was designed in [G87]. In the present paper we deal with a linear fragment of this general problem and describe

Typeset by AMS-TEX

for it an algorithm with a considerably better (than for the general problem) complexity, namely from the complexity class NC, i.e. with polynomial time and polylog depth (parallel time), moreover the algorithm produces a “triangular” basis for the space of solutions of the system.

A close problem to the one under consideration is solving linear system over the ring R of differential operators for which an algorithm was designed in [G91] (even for the case of the differential operators with the coefficients in many variables from ℚ(X_{1}, ⋯, X_{n})). The problem considered in the present paper is more subtle from the point of view of the allowed transformations of a system since for the sake of equivalence we may not multiply the equations of the system by the differential operators, as we could do in the case of the linear systems over R (see [G91] ).

Therefore, we need to carry out elementary transformations with the matrix over R of the system (see section 1 below), in order to reduce the matrix to a standard basis form which is a particular case of a differential standard basis [G], [O], [C] for partial differential operators. Since the ring R is non-commutative (some of its properties one can find in e.g. [B] ), the difficulties arise in estimating the standard basis form of the matrix over R unlike the case of the matrices over the (euclidean) rings of integers or univariate polynomials, because for the latter one exploits the notion of the determinant. But still we are able (see lemma 4) to bound the size of a quasiinverse of a matrix over R (for an inversible matrix a similar bound follows from the bound in [O] obtained for a more general situation of nonlinear operators) and define the rank of a matrix over R (see [J], also Lemma 5 below). To replace the notion of the determinant we consider (see section 2 below) the order [R] of a system of linear differential operators, i.e. of a matrix over R, being the dimension over ℂ(X) of the factor of the free R-module over the submodule generated by the rows of the matrix. We prove that the order is additive with respect to the product of the square matrices (Lemmas 6, 7). Relying on Lemma 7, on the analogue of Bezout’s theorem for differential equations [R], [Ko] (see also Lemma 9 below) and on a bound on a quasiinverse (see Lemma 4 in section 1), we estimate in section 3 the size of the standard basis form of the matrix (see Lemma 10) using the construction of a minimal element in a module with respect to a non-archemedian form (the order). In the last section 4 we give an algorithm from NC for constructing the standard basis form of a matrix, applying the bounds achieved in section 3. This provides a desired algorithm from NC for testing solvability of a system of linear ordinary differential equations and producing a “triangular” basis for the space of solutions of a system (see the theorem at the end of the paper).

Let us underline that the main purpose of this paper is to describe an algorithm with the low complexity (NC) for an important problem in symbolic computations in systems of linear differential equations. The needed auxiliary bounds from sections 1,2 (unfortunately, nowhere written explicitely) could be obtained without difficulties by the experts in differential algebra and they are included to make the paper self-contained.

Mention also that the problem of solving a single linear ordinary differential equation in one unknown leads to the problem of factoring of the equation, for the latter problem an algorithm was proposed in [G88]. A slight generalization of this problem is solving a first-order system of linear ordinary differential equations, an algorithm for reducing a matrix of this system to the block-triangular form was exhibited in [G90]. A connection of the first-order linear systems with the general linear systems considered in the present paper, is discussed in section 4 below.

1. \underline{Transformations and the rank of matrices over the ring of linear differ-} \underline{ential operators}.

Denote by D = \frac{d}{dX} , R = ℂ(X)[D] , and by F a Picard-Vessiot closure (see [K] ), i.e. any linear differential equation L = ( ∑_{0≤i≤n} f_{i}D^{i}) u = 0 with the coefficients f_{i} ∈ F and the leading coefficient ℓc(L) = f_{n} ≠ 0 has n linearly independent over ℂ solutions in F, and furthermore, a subfield of constants of F (i.e. the elements c ∈ F such that Dc = 0) coincides with ℂ.

We consider a problem of solvability in F of a system of linear ordinary differential equations in several unknowns

∑_{1≤j≤s} L_{ij}u_{j} = b_{i}, 1 ≤ i ≤ k  (1)

where L_{ij} ∈ ℚ(X)[D], b_{i} ∈ ℚ(X) and the solutions u_{1}, ⋯, u_{s} should be in F^{s}. For an operator L = ∑_{0≤i≤n} f_{i}D^{i} ∈ R with ℓc(L) = f_{n} ≠ 0 denote n = ord L and by deg L denote ∑_{0≤i≤n} deg_{X} f_{i}. Consider k × s matrix L = (L_{ij}), assume that ord L ≤ r, deg L ≤ d, deg (b_{i}) ≤ d, i.e. ord L_{ij} ≤ r, deg L_{ij} ≤ d for all i, j. Assume also that the bit-size of each (rational) coefficient of L_{ij}, b_{i} does not exceed M.

Consider now k×s matrix A = (A_{ij}) with the entries A_{ij} ∈ R, assume that ord (A_{ij}) ≤ r. As the ring R is left-euclidean, making elementary transformations over R with the rows, one can reduce A to the following standard basis form, see [J] (it is a particular case of a characteristic set [R] which is considered in [R] in nonlinear case, or of a differential standard basis [C], [G], [O] )

Q = ( \begin{matrix}0 & … & 0Q_{1p_{1}} \\ 0 & … & … & 0Q_{2p_{2}} \\ 0 & … & … & … & 0Q_{3p_{3}} & ∗ \\ ⋮ & & & & & ... \\ 0 & … & … & … & … & … & 0Q_{ℓ,p_{ℓ}} & … \\ & & & ◯ & & & ◯ \end{matrix} )  (2)

where p_{1} < p_{2} < ⋯ < p_{ℓ}, all the rows starting with (ℓ + 1)-th vanish . Let us admit also as an elementary transformation the multiplication (from the left) of a row by a nonzero element from ℂ(X). In other words, there is k × k matrix B = (B_{ij}) over R being a product of elementary matrices such that BA = Q. The rows of Q provide a triangular basis of a left R-module R^{k}A ⊂ R^{s} generated by the rows of the matrix A.

The next lemma and the corollary one can deduce from the results in [J].

Lemma 1. A square k × k matrix A over R is inversible from the left if and only if A equals to a product of elementary matrices.

Corollary. A square matrix is inversible from the left if and only if it is inversible from the right. The left inverse is unique and coincides with the right inverse. Thus, one could talk simply about inversible matrices.

We say that k ×s matrix A is quasiinversible from the left if there exists s×k matrix G over R such that GA = ( ..._{◯ C_{s}}^{C_{1} ◯} ) is a diagonal matrix with nonzero diagonal elements C_{1}, …, C_{s} (in a similar way one could define quasiinversibility from the right).

Lemma 2. A is quasiinversible from the left iff the dimension dim_{ℂ(X)}(R^{s}/R^{k}A) < ∞ of the factor-module is finite.

Proof. If GA = ( ..._{◯  C_{s}}^{C_{1}  ◯} ) and ord C_{1} = r_{1}, …, ord C_{s} = r_{s} then the vectors Π(e_{i}^{(j)}) =

Π( \underbrace{0, ⋯, 0, D^{j}}_{i}, 0, ⋯, 0) ∈ R^{s}/R^{k}A for 1 ≤ i ≤ s, 0 ≤ j < r_{i} constitute a generating set

over ℂ(X) of R-module R^{s}/R^{k}A, where Π : R^{s} → R^{s}/R^{k}A, is the natural projection, hence dim_{ℂ(X)}(R^{s}/R^{k}A) ≤ r_{1} + ⋯ + r_{s} (a better inequality see below in Lemmas 9, 10).

Let dim_{ℂ(X)}(R^{s}/R^{k}A) < ∞. Then one can reduce A by elementary transformations of the rows to standard basis form (2) and if ℓ < s then the infinite family of vectors Π(e_{p}^{(0)}), Π(e_{p}^{(1)}), …, where 1 ≤ p ≤ s is distinct from p_{1}, …, p_{ℓ}, are independent over ℂ(X) and we get a contradiction. Therefore, ℓ = s. One can show that there exists s × k matrix G over R such that GQ = ( ..._{◯ C_{s}}^{C_{1} ◯} ) with nonzero C_{1}, …, C_{s}. Indeed, multiply the first row of Q by a suitable element 0 ≠ α_{1} ∈ R such that α_{1}Q_{12} = β_{1}Q_{22} for a certain β_{1} ∈ R (this is possible since R is an Ore domain [B]), then subtract from the first row the second one multiplied by β_{1}, thereby we’ll achieve vanishing the entry with the coordinates (1, 2). Continuing in a similar way, we’ll make all the entries in the first row (except the diagonal entry) to be zeroes. Then we proceed to the second row and so on. As a result we’ll get a diagonal matrix which shows that A is quasiinversible from the left.

Observe that when dim_{ℂ(X)}(R^{s}/R^{k}A) < ∞, the latter dimension coincides with the order of the system Au = 0 [R]. In [R] the order was introduced for a prime ideal in the ring of differential polynomials, we use it for a linear ideal generated by the rows of A.

The next lemma was actually proved in [G91].

Lemma 3. A is quasiinversible from the left iff there does not exist a vector 0 ≠ v ∈ R^{s} such that Av = 0. For (s − 1) × s matrix A one can select 0 ≠ v ∈ R^{s} such that Av = O and ord (v) ≤ (s − 1)r + 1.

Proof. If A is quasiinversible from the left and GA = ( ..._{◯  C_{s}}^{C_{1}  ◯} ) then Av ≠ 0 for any

0 ≠ v ∈ R^{s} as R has no divisors of zero ( [B] ).

Conversely, let Av ≠ 0 for any 0 ≠ v ∈ R^{s}. Let us show that in the standard basis form (2) ℓ = s. Suppose ℓ < s. Consider the ℂ(X)-space R^{s;N} of the vectors (α_{1}, …, α_{s}) ∈ R^{s} for which ord (α_{1}), …, ord (α_{s}) < N. Let ord (Q_{ij}) ≤ R for all i, j. Then the composition of the mapping Q : v → Qv with the restriction onto first ℓ coordinates (notice that others are zeroes, see (2)) maps Q̄ : R^{s;N} → R^{ℓ;N+R}. As dim_{ℂ(X)} R^{s;N} = s N, for N = [\frac{ℓR}{s−ℓ}] + 1 we get that dim_{ℂ(X)} R^{s;N} > dim_{ℂ(X)} R^{ℓ;N+R} and therefore, there exists a vector 0 ≠ v ∈ R^{s;N} such that Qv = 0, hence BAv = 0 and thus Av = 0 since B is a product of elementary matrices (cf. Lemma 1). The obtained contradiction with the supposition justifies the equality ℓ = s. Then one can show that A is quasiinversible from the left. This proves the first statement of the lemma. For the second statement follow the latter proof considering instead of Q̄ the mapping Ā : R^{s;M} → R^{s−1;M+ord (A)} for M = (s − 1)ord A + 1.

The next lemma was proved in [G91].

Lemma 4. A square s × s matrix A is quasiinversible from the left iff A is quaisiinversible from the right. In this case there exists G for which GA = ( ..._{◯ C_{s}}^{C_{1} ◯} ) with ord (G) ≤ (s − 1)r + 1.

Proof. Let A be quasiinversible from the left. Then for an appropriate matrix B, being a product of elementary matrices, we have BA = ( ... ..._{◯ Q_{ss}}^{Q_{11} Q_{12}} ) where Q_{11} ⋯ Q_{ss} ≠ 0 (see (2) and the proof of Lemma 2). Let us show that for any vector 0 ≠ w ∈ R^{s} holds w A ≠ 0, this would imply that A is quasiinversible from the right because of Lemma 3. Assume that 0 = wA. Then 0 = wA = (wB^{−1}Q) and we get a contradiction, which justifies that A is quasiinversible from the right.

In order to prove the necessary bound on G, consider for each 1 ≤ j ≤ s a matrix A^{(j)} obtained from A by deleting its j-th column. Lemma 3 shows that there exists a vector 0 ≠ g^{(j)} ∈ R^{s} such that g^{(j)}A^{(j)} = 0 and ord g^{(j)} ≤ (s − 1)r + 1. As a matrix G take a matrix with j-th row equal to g^{(j)}.

Notice that when A is inversible, lemma 4 follows from the theorem 6 in [O], where a similar bound was proved for a much more general situation of an inversible nonlinear differential map.

Thus, for a square matrix A we can say that it is quasiinversible without specifying from the left or from the right. Notice (see also [G91] ) that a square matrix A is quasiinversible iff its Dieudonné determinant ( [A] ) does not vanish.

Define the rank rk(A) as a maximal ℓ such that there exists ℓ × ℓ quasiinversible submatrix of A (cf. [J] ), the following lemma can be deduced from the results in [J].

Lemma 5. rk(A) coincides with a) ℓ in the standard basis form (2);

b) the maximal number of the columns of A being R-linearly independent;

c) the maximal number of the rows of A being R-linearly independent.

2. \underline{Some properties of the order of a system of linear differential operators}. For brevity we adopt the notation dim(R^{s}/R^{k}A) = dim_{ℂ(X)}(R^{s}/R^{k}A).

Lemma 6. For any m × k matrix B and k × s matrix A over R dim(R^{s}/R^{m}BA) ≤ dim(R^{k}/R^{m}B) + dim(R^{s}/R^{k}A)

If A is quasiinversible from the right then this inequality turns to be the equality.

Proof. Consider the natural projections Π_{1} : R^{s} → R^{s}/R^{k}A, Π_{2} : R^{k} → R^{k}/R^{m}B, Π_{3} : R^{s} → R^{s}/R^{m}BA. Let u_{1}, …, u_{γ} ∈ R^{k} be such that Π_{2}(u_{1}), …, Π_{2}(u_{γ}) constitute a basis over ℂ(X) of R^{k}/R^{m}B, and v_{1}, …, v_{δ} ∈ R^{s} be such that Π_{1}(v_{1}), …, Π_{1}(v_{δ}) constitute a basis over ℂ(X) of R^{s}/R^{k}A (note that γ or δ could be infinite). Let us prove that Π_{3}(v_{1}), …, Π_{3}(v_{δ}), Π_{3}(u_{1}A), …, Π_{3}(u_{γ}A) generate R^{s}/R^{m}BA over ℂ(X) and constitute a basis when A is quasiinversible from the right. Indeed, let for some elements f_{1}, …, f_{δ}, g_{1}, …, g_{γ} ∈ ℂ(X) and a vector (β_{1}, …, β_{m}) ∈ R^{m} we have f_{1}v_{1} + ⋯ + f_{δ}v_{δ} + (g_{1}u_{1} + ⋯ + g_{γ}u_{γ})A = (β_{1}, ⋯, β_{m})BA, then f_{1} = ⋯ = f_{δ} = 0. If A is quasiinversible from the right then g_{1}u_{1} + ⋯ + g_{γ}u_{γ} = (β_{1}, ⋯, β_{m})B by virtue of Lemma 3, hence g_{1} = ⋯ = g_{γ} = 0.

On the other hand, for any vector w ∈ R^{s} there exist f_{1}, …, f_{δ} ∈ ℂ(X) and a vector v ∈ R^{k} for which w = f_{1}v_{1} + ⋯ + f_{δ}v_{δ} + vA. Then v = g_{1}u_{1} + ⋯ + g_{γ}u_{γ} + uB for suitable g_{1}, …, g_{γ} ∈ ℂ(X), u ∈ R^{m}. Therefore w = f_{1}v_{1} +⋯+f_{δ}v_{δ} +g_{1}u_{1}A+⋯+g_{γ}u_{γ}A+uBA, i.e. dim(R^{s}/R^{m}BA) ≤ γ + δ = dim(R^{k}/R^{m}B) + dim(R^{s}/R^{k}A).

In other terms we can reformulate what was proved above, saying that we have the following exact sequence of ℂ(X)-vector spaces

R^{k}/R^{m}B ⟶^{α} R^{s}/R^{m}BA ⟶^{π} R^{s}/R^{k}A → O

where α(v + R^{m}B) = vA + R^{m}BA and π(w + R^{m}BA) = w + R^{k}A. In the case of quasiinversible A the following sequence is exact:

O → R^{k}/R^{m}B ⟶^{α} R^{s}/R^{m}BA ⟶^{π} R^{s}/R^{k}A → O

Lemma 7. If a matrix A is square then dim(R^{s}/R^{m}BA) = dim(R^{s}/R^{m}B)+dim(R^{s}/R^{s}A)

Proof. If A is quasiinversible (see Lemma 4) then we use Lemma 6.

If A is not quasiinversible then dim(R^{s}/R^{m}BA) ≥ dim(R^{s}/R^{s}A) = ∞.

Remark that as in the following example dim (R^{2}/R^{2} ( \begin{matrix}1 & 1 \\ 1 & D \end{matrix} )) = 1, dim (R/R^{2} ( \begin{matrix}1 \\ 1 \end{matrix} )) = 0 and for the product of these matrices dim (R/R^{2}\binom{2}{1+D}) = 0, the inequality in Lemma 6 for rectangular matrices could be strict.

Lemma 8. a) For a triangular k × s (where k ≥ s) matrix C = ( ..._{◯  C_{s}}^{C_{1}  ∗} ) we have

dim(R^{s}/R^{k}C) = ord C_{1} + ⋯ + ord C_{s}, provided that C_{1} ⋯ C_{s} ≠ 0.

b) dim(R^{s}/R^{k}A) < ∞ iff ℓ = s in the standard basis form (2). In this case dim(R^{s}/R^{k}A) = ord Q_{11} + ⋯ + ord Q_{ss}.

c) When A is a square matrix then dim(R^{s}/R^{s}A) = dim(R^{s}/AR^{s}), where in the right side of the equality we regard R^{s} as a right R-module.

d) A square matrix A is inversible iff dim(R^{s}/R^{s}A) = 0.

Proof. a) Is obvious.

b) The first statement one can find in the proof of Lemma 3. The second statement follows from a) and the equality dim(R^{s}/R^{k}A) = dim(R^{s}/R^{k}Q).

c) Because of Lemma 4 both left and right sides of the equality are finite or infinite simultaneously. Assume they are both finite. Then BA = ( ..._{◯ Q_{ss}}^{Q_{11} ∗} ) (see (2)) where B is a product of elementary matrices and Q_{11} ⋯ Q_{ss} ≠ 0 (see b) ). For any s×s elementary matrix G we have dim(R^{s}/R^{s}G) = dim(R^{s}/GR^{s}) = 0, hence by Lemma 7 the same is true for any inversible matrix (cf. Lemma 1), thus dim(R^{s}/R^{s}B) = dim(R^{s}/BR^{s}) = 0.

a) implies that for the triangular matrix

Q = ( ..._{◯ Q_{ss}}^{Q_{11} ∗} ) the equalities dim(R^{s}/R^{s}Q) = dim(R^{s}/QR^{s}) = ord Q_{11} + ⋯ +

ord Q_{ss} hold, then Lemma 7 entails c).

d) follows from b) and Lemma 1.

The following lemma was proved in [R], p. 135 (see also [Ko] ) in a more general form for the order of a prime ideal in the ring of differential polynomials.

Lemma 9. If k × s matrix A is quasiinversible from the left then dim(R^{s}/R^{k}A) ≤

max_{i} {ord a_{i1}} + ⋯ + max_{i} {ord a_{is}}.

3. \underline{Bounds on the standard basis form of a matrix over the ring of differen-} \underline{tial operators}.

In this section we’ll estimate ord (Q), ord (B) in the standard basis form (2) relying on the results on the order from the section 2.

Take any s×s permutation matrix P, mapping P (p_{1}) = 1, …, P (p_{ℓ}) = ℓ, then BAP =

( \begin{matrix}Q_{1p_{1}} \\ & . & . & . \\ & & & & Q_{ℓp_{ℓ}} & ... \\ ◯ & & & & & ◯ \end{matrix} ). Represent AP = (A_{1}A_{2}) where k × ℓ submatrix A_{1} consists of the

first ℓ columns of AP, then by Lemma 5 rkA = rkA_{1} = ℓ. Complete A_{1} by (k −ℓ) columns of the type (0, …, 0, 1, 0, …, 0)^{T} to k × k quasiinversible matrix (A_{1}A_{3}). Then

B(A_{1}A_{3}) = ( \begin{matrix}Q_{1p_{1}} \\ . & . & . \\ & & Q_{ℓp_{ℓ}} \\ & ◯ & & ∗ \end{matrix} )

Making several elementary transformations with the rows having indices bigger than ℓ, reduce the matrix at the right side to the triangular form

B_{0}(A_{1}A_{3}) = ( \begin{matrix}Q_{1p_{1}} & & & ∗ \\ . & . & . \\ & & Q_{ℓp_{ℓ}} \\ & & & Q_{ℓ+1,ℓ+1}^{(0)} \\ & ◯ & & ... \\ & & & & Q_{kk}^{(0)} \end{matrix} )

herewith B = ( \begin{matrix}B_{1} \\ B_{2} \end{matrix} ) where B_{1} is ℓ × k submatrix, B_{0} = ( \begin{matrix}B_{1} \\ B_{3} \end{matrix} ) and dim(R^{k}/R^{k}B) = dim(R^{k}/R^{k}B_{0}) = 0 (see Lemma 8 d)).

Moreover, making some elementary transformations with the rows, one can assume w.l.o.g. that ord (Q_{ip_{j}}) < ord (Q_{jp_{j}}), ord (Q_{ij}^{(0)}) < ord (Q_{jj}^{(0)}) for all i < j.

By Lemmas 6, 7, 9 ord (Q_{1p_{1}}) + ⋯ ord (Q_{ℓp_{ℓ}}) + ord (Q_{ℓ+1,ℓ+1}^{(0)}) + ⋯ + ord (Q_{kk}^{(0)})

= dim(R^{k}/R^{k}A_{1}A_{3}) ≤ max_{i}{ord A_{ip_{1}}} + ⋯ + max_{i}{ord A_{ip_{ℓ}}} ≤ ℓr, hence ord (Q_{i,p_{j}}),

ord (Q_{ij}^{(0)}) ≤ ℓr. By Lemma 4 there exists k × k matrix G over R such that (A_{1}A_{3}) G =

( ..._{◯  C_{k}}^{C_{1}  ◯} )
where C_{1} ⋯ C_{k} ≠ 0 and ord (G) ≤ (k−1)r+1, hence ord (C_{i}) ≤ kr+1. As B_{0} ( ..._{◯  C_{k}}^{C_{1}  ◯} )

= ( \begin{matrix}Q_{1p_{1}} & & & & & & & ∗ \\ & . & . & . \\ & & & & . & . & . \\ ◯ & & & & & & & Q_{kk}^{(0)} \end{matrix} ) G, we conclude that ord (B_{0}) ≤ (ℓ + k − 1)r + 1.

Observe that B_{0}A has the standard basis form similar to (2) (with the same “diagonal” entries Q_{1p_{1}}, …, Q_{ℓ,p_{ℓ}} and perhaps different other entries as we achieved the conditions ord Q_{ip_{j}} < ord Q_{jp_{j}}, ord (Q_{ij}^{(0)}) < ord (Q_{jj}^{(0)}), i < j)

B_{0}A = ( \begin{matrix}0 & … & 0 Q_{1p_{1}} & & & ∗ \\ 0 & … & … & 0 Q_{2p_{2}} \\ ⋮ & & & & ... \\ 0 & … & … & … & … & 0 Q_{ℓp_{ℓ}} \\ & & & & ◯ & & ◯ \end{matrix} )

since rkA = ℓ. Therefore ord (Q) ≤ (ℓ + k)r + 1. Let us summarize the proved in the present section in the following lemma.

Lemma 10. There exists an inversible matrix B such that BA = Q has the standard basis form (2) and moreover ord (B) ≤ (s + k − 1)r + 1, ord (Q) ≤ (s + k)r + 1.

4. \underline{NC algorithm for finding standard basis form of a matrix over the ring} \underline{of differential operators}.

Let us design an algorithm which finds the standard basis form of a matrix in NC, i.e. polynomial time and with polylogarithmic depth (parallel complexity).

Join to the matrix A the unit matrix and denote the resulting k × (s + k) matrix by Ā = (AE). Obviously rkĀ = k (see Lemma 5). Therefore, the standard basis form of Ā equals to

B_{1}Ā = ( \begin{matrix}0 & … & 0 Q_{1p_{1}} \\ ⋮ & & & ... \\ 0 & … & … & … & 0 Q_{ℓp_{ℓ}} & ∗ \\ ⋮ & & & & & ... \\ 0 & … & … & … & … & … & 0 Q_{kp_{k}} & … \end{matrix} ) = Q̄ (see (2))

where dim(R^{k}/R^{k}B_{1}) = 0 (see Lemma 8).

For each 1 ≤ m ≤ s + k and 0 ≤ j ≤ (s + 2k)r + 1 the algorithm tests, whether there exists a vector w = (w_{1}, …, w_{k}) with ord (w) ≤ (s+2k −1)r +1 (cf. Lemma 10) such that the vector wĀ = ( \underbrace{0, …, 0, v}_{m}, …) where ord v = j and the leading coefficient ℓc(v) = 1. The latter condition can be written as a linear system T_{m,j} over ℚ(X) with k((s+2k−1)r+ 2) unknowns being the coefficients of w_{1}, …, w_{k} in the powers of 1, D, …, D^{(s+2k−1)r+1} and with at most (s + k)((s + 2k)r + 1) equations. As the entries of this linear system are the rational functions from ℚ(X) with the degrees in X not exceeding d and with the size of rational coefficients at most M, the algorithm can solve T_{m,j} in time (Mdskr)^{0(1)} with the depth (parallel complexity) log^{0(1)}(Mdskr) using [M].

As the rows of the matrix Q̄ constitute a (triangular) basis of the left R-module R^{k}Ā, the system T_{m,j} could be solvable only for m = p_{1}, …, p_{k}. For each of these m = p_{i} take the minimal such j_{i} for which T_{p_{i},j_{i}} is solvable. Lemma 10 implies that T_{p_{i},j} is solvable for j = ord Q_{i,p_{i}}, hence j_{i} ≤ ord Q_{i,p_{i}}. Take a solution W^{(i)} = (w_{1}^{(i)}, …, w_{k}^{(i)}) for T_{p_{i},j_{i}}, and denote by W k × k matrix with i-th row to be W^{(i)}. Then

W Ā = ( \begin{matrix}0 & … & 0 Q̃_{1p_{1}} \\ ⋮ & & & ∗ \\ 0 & … & … & 0 Q̃_{kp_{k}} & … \end{matrix} ) = Q̃

where ord Q̃_{i,p_{i}} = j_{i}.

Let us prove that dim(R^{k}/R^{k}W) = 0. Denote by Ā_{0}, Q̄_{0}, Q̃_{0} the k × k matrices obtained from Ā, Q̄, Q̃ respectively, by taking the submatrices formed by the columns p_{1}, …, p_{k}. Then B_{1}Ā_{0} = Q̄_{0}, W Ā_{0} = Q̃_{0}.

Lemmas 7, 8 a) entail 0 = dim(R^{k}/R^{k}B_{1}) = ord Q_{1,p_{1}}+⋯+ord Q_{k,p_{k}}− dim(R^{s+k}/R^{k}Ā) ≥ ord Q̃_{1,p_{1}}+⋯+ord Q̃_{k,p_{k}}− dim(R^{s+k}/R^{k}Ā) = dim(R^{k}/R^{k}W) ≥ 0, therefore dim(R^{k}/R^{k}W) = 0 and moreover ord Q̃_{i,p_{i}} = ord Q_{i,p_{i}}, 1 ≤ i ≤ k. As W A has a desired standard basis form (see (2)), we get the following lemma

Lemma 11. There is an NC-algorithm, so running in time (Md skr)^{0(1)} with a depth (parallel complexity) log^{0(1)}(Md skr), which produces an inversible over R k × k matrix

W such that W A = ( \begin{matrix}0 & ... & 0 Q̃_{1p_{1}} \\ ⋮ & & & ∗ \\ 0 & ... & ... & 0 Q̃_{ℓp_{ℓ}} \\ & & ◯ & & ◯ \end{matrix} ) has the standard basis form.

Now we get a criterium for solvability of a system (1). Namely, apply Lemma 11 to k × (s + 1) matrix A = (L_{ij} b_{i})_{1≤i≤k, 1≤j≤s}, so the last column is ( ⋮_{b_{k}}^{b_{1}} ). Then the system

(1) has a solution in the field F iff and only if p_{ℓ} ≤ s (in other words p_{ℓ} ≠ s + 1) and the standard basis form provides a “triangular” basis of the space of solutions of (1). Let us summarize the obtained above in the following main result of the paper.

Theorem. One can test solvability of a system (1) of linear differential equations in several unknowns in the Picard-Vessiot closure F and find a “triangular” basis of the space of solutions of (1) in NC complexity class, so with the time (Md skr)^{0(1)} and with a depth (parallel time) log^{0(1)}(Md skr).

Observe that the space of solutions of a homogeneous system (1), so when b_{1} = ⋯ = b_{k} = 0, has a finite dimension (over ℂ(X)) if and only if p_{1} = 1, …, p_{ℓ} = ℓ and ℓ = s (for k × s matrix A = (L_{ij}), see above). In this case the standard basis form W A of the system can be rewritten in the common first-order matrix form DY = HY (cf. [G90] ) where the vector Y has coordinates u_{1}, Du_{1}, …, D^{j_{1}−1}u_{1}, u_{2}, …, D^{j_{2}−1}u_{2}, …, u_{s}, …, D^{j_{s}−1}u_{s} and j_{i} = ord Q̃_{i,p_{i}}, 1 ≤ i ≤ s, one could easily get the matrix H over ℚ(X) from the matrix W A.

Acknowledgements. The author would like to thank Mike Singer for the attention to the paper.

REFERENCES

[A] Artin, E., Geometric algebra, Interscience Publ., 1957.

[B] Björk, J.-E., Rings of differential operators, North-Holland, 1979.

[C] Carro-Ferro, G., Groebner bases and differential ideals, Lect. Notes Comput. Sci., 356 (1987), pp. 129–140.

[G] Galligo, A., Some algorithmic questions on ideals of differential operators, Lect. Notes Comput.

Sci., 204 (1985), pp. 413–421.

[G87] Grigoriev, D., Complexity of quantifier elimination in the theory of ordinary differential equations, Lect. Notes Comput. Sci., 378 (1989), pp. 11–25.

[G88] Grigoriev, D., Complexity of factoring and GCD calculating of ordinary linear differential operators, J. Symb. Comput., 10, N1 (1990), pp. 7–37.

[G90] Grigoriev, D., Complexity of irreducibility testing for a system of linear ordinary differential equations, Proc. Int. Symp. on Symb. Algebr. Comput., ACM (1990) Japan, pp. 225–230.

[G91] Grigoriev, D., Complexity of solving systems of linear equations over the rings of differential operators, Proc. Int. Symp. Eff. Meth. in Algebraic Geometry (1990) Italy; in Progr. in Math., Birkhäuser, v.94, 1991, pp. 195–202 .

[J] Jacobson, N., Pseudo-linear transformations, Ann. Math., 38, N2 (1937), pp. 484–507.

[K] Kaplansky, I., An introduction to differential algebra, Hermann, Paris, 1957.

[Ko] Kolchin, E., Differential algebra and algebraic groups, Academic Press, 1973.

[M] Mulmuley, K., A fast parallel algorithm to compute the rank of a matrix over an arbitrary field, Proc. 18 STOC ACM (1986), pp. 338–339.

[O] Ollivier, F., Standard bases of differential ideals, Lect. Notes Comput. Sci., 508 (1991), pp. 304–321.

[R] Ritt, J.F., Differential algebra, Amer. Math. Soc. Colloq. Publ., vol. 33, NY, 1950.

Departments of Computer Science and Mathematics Pennsylvania State University University Park, PA 16802 USA e-mail: dima@math.psu.edu