       Re: Calculating the null space of a matrix with univariate polynomial entries?

• To: mathgroup at smc.vnet.net
• Subject: [mg67260] Re: [mg67227] Calculating the null space of a matrix with univariate polynomial entries?
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Wed, 14 Jun 2006 06:29:48 -0400 (EDT)
• References: <200606130507.BAA23813@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

```scott.morrison at gmail.com wrote:
> Hi,
>
> I'm having great difficulty calculating null spaces of (fairly sparse)
> matrices with univariate polynomial entries. There's an example below.
>
> First, I'm finding NullSpace[m, Method->OneStepRowReduction] gives much
> better results than any of the other methods, but it still very
> unsatisfactory. My typical experience with larger matrices (roughly
> 40x40), such as the one below, is that Mathematica will work for quite
> some time (25 minutes, in the example below), and then start allocating
> huge amounts of memory (hundreds of MBs), causing paging and an
> eventual crash.
>
> This doesn't seem reasonable, although I haven't yet tried these
> examples in any other computer algebra systems.
>
> Does anyone know what's going on? Is there something I can do to avoid
> this issue? Does anyone know of an implementation of NullSpace specific
> to univariate polynomial matrices?
>
> Any help much appreciated!
>
> Scott Morrison,
> UC Berkeley
>
> This example is also online at
> http://math.berkeley.edu/~scott/A1_2_5_0.m in case it gets munged here.
> Incidentally, I know the null space is meant to be 6 dimensional, if
> that happens to help (and I'll always know this dimension for my other
> examples).
> [matrix removed for brevity --dl]

This comes up from time to time. My recommendation is to do it by
interpolation. I would go about it as follows. Clear denominators,
transpose, tack on an identity matrix to the right (in essence to record
relations), and row reduce with q set to various integer values. After
that one must figure out which reductions gave a null space that will be
correct on interpolating (that is, which ones do not have extra null
rows), and use the corresponding evaluations for interpolation.

The code below is not bullet proof. I do not test a priori that I have
sufficiently many interpolating values. One could base such a number on
an estimate of maximal degrees of minors. Also I tacitly assume that the
interpolation should give a polynomial as opposed to a rational
function. Maybe this is justifiable, I'm not sure offhand, though
certainly by clearing denominators one can obtain null vectors that are
polynomial.

(* Step 1: Clear denominators, transpose, append an identity matrix. *)

m2 = Together[mat];
denom = Apply[PolynomialLCM,Flatten[Denominator[m2]]]
m3 = m2*denom;
len = Length[m3[]];
l2 = Length[m3];
m4 = Transpose[Join[m3, IdentityMatrix[len]]];

(* Step 2: Do row reductions for evaluated matrix. *)

rforms = Table[RowReduce[m4], {q,-10,10}];

(* Step 3: Extract rows corresponding to null vectors. Find the correct
number of null vectors. *)

hasNull[vec_,n_] := Take[vec,n]===Table[0,{n}]
nulls = Map[Select[#,hasNull[#,l2]&]&, rforms];
numnulls = Min[Map[Length,nulls]]

(* Step 4: Figure out which reductions gave the right number of null
vectors, and which are the corresponding  interpolation values. *)

order = Ordering[Map[Length,nulls]];
nulls2 = Select[nulls, Length[#]===numnulls&];
interpvals = Take[Range[-10,10][[order]], Length[nulls2]]

(* Step 5: Remove row prefixes corresponding to zeroed elements.
Transpose to put into convenient form for interpolation. *)

nulls3 = Map[Drop[#,l2]&, nulls2, {2}];
nulls4 = Transpose[nulls3,{3,1,2}];

(* Step 6: Interpolate. *)

interp[vec_] := InterpolatingPolynomial[Transpose[{interpvals,vec}],q]
nulls5 = Map[Expand[interp[#]]&, nulls4, {2}];

Now check the result. First see the dimensions are correct.

In:= Dimensions[nulls5]
Out= {6, 51}

Here are the null vectors.

In:= InputForm[nulls5]
Out//InputForm=

{{1, 0, -q^4, 0, 0, -q^4, 0, 0, q^6 + q^8, 0, 0, -q^8, q^10, q^6 - q^10,
q^8,
-q^8 - q^10, -q^4, -q^4, -q^2 + q^6 + q^8, 0, q^6, q^4 - q^8 - q^10, 0,
q^6, q^10, q^6 - q^12, q^4 - 2*q^8 - q^10 + q^12, -q^10,
-q^6 - q^8 + q^10 + q^12, q^8, -q^10, -q^6 + q^10, -q^8, -q^6 + q^8 +
q^10,
q^6, q^6 + q^8, q^4 - q^8 - q^10, q^6 - q^10, q^4 - 2*q^8 - q^10 + q^12,
q^2 - q^4 - 2*q^6 + q^8 + 2*q^10, -q^6 + q^10, -q^4 + 2*q^8 - q^12,
-q^8 - q^10, -q^6 - q^8 + q^10 + q^12, -q^4 + 2*q^8 - q^12,
-q^6 + q^8 + q^10, -q^4 + q^6 + 2*q^8 - q^10 - q^12, q^4 - q^8, q^6,
q^4 - q^8, -q^6}, {0, 1, -q^2, 0, 0, 0, 0, -q^4, q^4 + q^6, 0, 0, -q^6,
q^8, -q^8, q^6, -q^8, 0, -q^4, q^6, -q^2, -q^2 + q^4 + q^6,
q^2 - q^6 - q^8, 0, 0, q^8, q^6 - q^10, -q^6 - q^8 + q^10, q^4 - q^8,
-q^6 + q^10, q^6, -q^8, q^8, -q^6, q^8, 0, q^6, -q^8, q^4 + q^6 - q^8,
q^4 - q^6 - 2*q^8 + q^10, -q^4 + 2*q^8, q^2 - q^4 - q^6 + q^8,
-q^4 + q^6 + q^8 - q^10, -q^6 - q^8, -q^6 + q^8 + q^10, q^6 - q^10,
-q^4 + q^8, q^6 - q^10, 0, q^4, -q^6, 0}, {0, 0, 0, 1, 0, -q^4, 0, 0,
q^6,
0, -q^2, -q^6, q^8, q^4 + q^6 - q^8, q^6, -q^6 - q^8, -q^2, 0, q^6, 0, 0,
-q^8, -q^4, -q^2 + q^4 + q^6, q^8, q^6 - q^10, q^4 - q^6 - 2*q^8 + q^10,
-q^8, -q^6 + q^8 + q^10, q^6, q^4 - q^8, q^2 - q^4 - q^6 + q^8, -q^6,
-q^4 + q^8, q^4, q^4 + q^6, q^2 - q^6 - q^8, -q^8, -q^6 - q^8 + q^10,
-q^4 + 2*q^8, q^8, q^6 - q^10, -q^8, -q^6 + q^10, -q^4 + q^6 + q^8 -
q^10,
q^8, q^6 - q^10, -q^6, 0, 0, 0}, {0, 0, 0, 0, 1, -q^2, 0, -q^2, q^2 +
q^4,
0, 0, -q^2 - q^4, -q^2 + q^4 + q^6, q^2 - q^6, q^2 + q^4, -q^4 - q^6, 0,
-q^2, q^4, 0, q^4, -q^4 - q^6, -q^2, q^4, -q^2 + q^4 + q^6,
-q^2 + 3*q^4 - q^6 - q^8, q^2 - 2*q^4 - q^6 + q^8, q^2 - q^4 - q^6,
-q^4 + q^6 + q^8, q^2 + q^4, q^2 - q^4 - q^6, -q^2 + q^6, -q^4, q^6, 0,
q^2 + q^4, -q^4 - q^6, q^2 - q^6, q^2 - 2*q^4 - q^6 + q^8,
-q^2 + q^4 + 2*q^6, -q^2 + q^6, q^4 - q^8, -q^4 - q^6, -q^4 + q^6 + q^8,
q^4 - q^8, q^6, -q^8, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 1, -q^2, q^2, 0, 0,
-q^2, q^4, -q^4, 0, 0, 0, 0, 0, -q^2, q^4, -q^4, 0, 0, -q^2 + q^4,
q^4 - q^6, -q^4 + q^6, 0, 0, q^2, -q^4, q^4, 0, 0, 0, 0, 0, q^2,
-q^4, q^4,
0, 0, -q^4, q^6, -q^6, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
-q^2, -q^2, -q^2 + q^4, q^2, q^2, -q^4, 0, 0, 0, 0, 0, 0, -q^2, q^4, q^4,
q^4 - q^6, -q^4, -q^4, q^6, 0, 0, 0, 0, 0, 0, q^2, -q^4, -q^4, -q^4 +
q^6,
q^4, q^4, -q^6, 0, 0, 0, 0, 0, 0, 0, 0, 0}}

Now observe we get zeros from appropriate dot products.

In:= Select[Flatten[Map[Expand[m3.#]&, nulls5]], #=!=0&]
Out= {}

Daniel Lichtblau
Wolfram Research

```

• Prev by Date: Re: Determining continuity of regions/curves from inequalities
• Next by Date: Re: structure array equivalent in Mathematica
• Previous by thread: Calculating the null space of a matrix with univariate polynomial entries?
• Next by thread: scaled Bessel functions in Mathematica?