Re: A few questions about RowReduce - bis
- To: mathgroup at smc.vnet.net
- Subject: [mg24256] Re: [mg24074] A few questions about RowReduce - bis
- From: Murray Eisenberg <murray at math.umass.edu>
- Date: Wed, 5 Jul 2000 23:10:37 -0400 (EDT)
- Organization: Mathematics & Statistics, Univ. of Mass./Amherst
- References: <8j9dsq$51v@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
A reason to find a row-echelon form that is not the canonical "reduced"
row-echelon form, is that you can answer questions about the nature of
solutions of linear systems, or linear independence of vectors, etc.,
without having to do -- or asking Mathematica to do -- all the extra
work to get the reduced row-echelon form.
Preston Nichols wrote:
>
> Jack-Michel,
>
> Here are some ideas which I hope will be helpful, though they do not
> address your questions about a "non-reduced" or "partial" row-echelon form.
> (To be honest, I can't think of why that specific form would be
> essential.) My discussion focuses on equations for the image of a matrix.
>
> Here's your matrix:
>
> In[1]:= n = 5;
> m = Table[(i + j - 1)^2, {i, n}, {j, n}];
>
> A vector of variables (note I changed yi from you message to y[i]):
>
> In[2]:= v = Table[y[i], {i, n}];
>
> The same structure can be built this way:
>
> In[3]:= v = Array[y,n];
>
> The vector v appended to m as the last column:
>
> In[4]:= << "LinearAlgebra`MatrixManipulation`";
> mv = AppendRows[m, Transpose[{v}]]
>
> The same operation, without requiring the MatrixManipulation package:
>
> In[5]:= mv = Transpose[Append[Transpose[m], v]]
>
> RowReduce[mv] permits division by polynomials in the y-variables, even
> though such expression may have the value zero. Row-reduction of mv,
> treating the y-variables as generic, leads to (from your own Mathematica
> function):
>
> {{1, 4, 9, 16, 25, y[1]},
> {0, -7, -20, -39, -64, -4y[1] + y[2]},
> {0, 0, 8/7, 24/7, 48/7, (17*y[1])/7 - (20*y[2])/7 + y[3]},
> {0, 0, 0, 0, 0, -y[1] + 3*y[2] - 3*y[3] + y[4]},
> {0, 0, 0, 0, 0, -3*y[1] + 8*y[2] - 6*y[3] + y[5]}}
>
> from which can be extracted equations for the image of the matrix:
>
> -y[1] + 3 y[2] - 3 y[3] + y[4]==0,
> -3 y[1] + 8 y[2] - 6 y[3] + y[5]==0
>
> One of your questions is how to get Mathematica to perform such a
> calculation. I found three ways, each having its own interest.
>
> ====================================================
>
> Mathematica solution 1: The "Left Nullspace"
>
> In[6]:= leftnullspace = NullSpace[Transpose[m]]
>
> Out[6]= {{-3, 8, -6, 0, 1}, {-1, 3, -3, 1, 0}}
>
> These vectors span the orthogonal complement of the image of m, so by
> setting dot products with v to zero, you can recover the desired equations.
>
> In[7]:= Thread[leftnullspace.v == 0]
>
> Out[7]= {-3 y[1] + 8 y[2] - 6 y[3] + y[5] == 0,
> -y[1] + 3 y[2] - 3 y[3] + y[4] == 0}
>
> (Thread changes an "equation of lists" into a "list of equations".)
>
> This strikes me as "the right way", because (1) it is brief, using
> high-level built-in Mathematica functions, and (2) it relies on a
> fundamental theorem of finite-dimensional linear algebra. An aphorism:
> Sometimes Mathematica does mathematics for you; sometimes you must do
> mathematics for Mathematica.
>
> ====================================================
>
> Mathematica solution 2: Use Eliminate
>
> In[8]:= ImageEquations[mat_?MatrixQ, y_Symbol] := Module[{m, n, x},
> {m, n} = Dimensions[mat];
> Eliminate[Thread[mat.Array[x, n] == Array[y, m]], Array[x, n]]
> ]
>
> In[9]:= ImageEquations[m, y]
>
> Out[9]= y[1] == 3 y[2] - 3 y[3] + y[4] && y[2] == 3 y[3] - 3 y[4] + y[5]
>
> Note that this method returns equations, rather than vectors (of
> coefficients for the equations) as my other solutions do.
>
> To understand how this function works, examine the step-by-step output from
> input lines such as
>
> u = Array[x, n]
> v = Array[y, n]
>
> m.u == v
>
> Thread[%] (* "equation of lists" to "list of equations" again *)
>
> Eliminate[%, u]
>
> This may not always give exactly the same equations as the "left nullspace"
> method, but the difference is due only to a change of basis for the
> orthogonal space to the image of m. Mathematically, the answers are
> equivalent.
>
> ====================================================
>
> Mathematica solution 3: Use RowReduce after all!
>
> Here's a function which appends an IdentityMatrix of appropriate size to
> the right of the given matrix, and then applies RowReduce to the whole
> thing. (This is what one does when using row operations to find the
> inverse of a square matrix, but it can be done, in principle, on any matrix.)
>
> Here's a version that makes use of LinearAlgebra`MatrixManipulation`:
>
> In[10]:= doubleRowReduce[mat_?MatrixQ] :=
> RowReduce[
> AppendRows[mat, IdentityMatrix[First[Dimensions[mat]]]]
> ]
>
> Here's a version of the same function that doesn't require the package:
>
> In[11]:= doubleRowReduce[mat_?MatrixQ] :=
> RowReduce[
> MapThread[Join, {mat, IdentityMatrix[First[Dimensions[mat]]]}]
> ]
>
> The coefficients for the equations of the image of the matrix m appear at
> the right in the bottom rows of doubleRowReduce[m]. The rows to use are
> those which reduce to zeroes in m. (I won't try to reproduce the output
> here; try it and see!) You can think of the columns of the appended
> identity matrix as holding coefficients of the y-variables. The original
> problem required a cummulative "record" of the arithmetical operations
> performed on each of the y[i]'s, so giving each y[i] a column of its own
> makes sense.
>
> Here's a function to extract the desired coefficients from the results of
> doubleRowReduce, using TakeMatrix from LinearAlgebra`MatrixManipulation. I
> suspect this could be done in a better way, but it works:
>
> In[12]:= imagecoefficients[mat_?MatrixQ] :=
> With[{m = Length[mat], n = Length[First[mat]]},
> TakeMatrix[mat,
> {-Count[TakeColumns[mat, n - m], _?(Union[#] == {0} &)],
> n - m + 1},
> {-1, -1}
> ]
> ]
>
> In[13]:= imagecoefficients[doubleRowReduce[m]]
>
> Out[13]= {{1, 0, -6, 8, -3}, {0, 1, -3, 3, -1}}
>
> Again, this may not always give exactly the same equations as the other
> methods, but the difference is due only to a change of basis in the
> y-variables, leaving the image and its complementary space invariant.
>
> ====================================================
>
> I hope this is of some use for you.
>
> Preston Nichols
> Mathematics and Computer Science
> Wittenberg University
>
> At 02:26 AM 06/23/2000 -400, jmcornil wrote:
> >First question : >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
> >
> >The function RowReduce perform a Gauss-Jordan reduction putting
> >1 on the "diagonal" and 0 elsewhere so long it is possible. But is there
> >an other function (or an option for RowReduce) performing a simple
> >Gauss reduction i.e. putting 0 ONLY UNDER the first diagonal ?
> >
> >By the way I saw there is an option "Method" for the function.
> >
> >In[22]:= Options[RowReduce]
> >Out[22]= {Method -> Automatic, Modulus -> 0, ZeroTest -> Automatic}
> >
> >and I did not found any help on it. How to use it ?
> >
> >Second question : >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
> >
> >I did not foud any function to perform a "limited" Gauss reduction, i.e.
> >a function performing
> >a reduction up to a column number . Such an option woud be useful in a
> >case as follows
> >to get equations of the image of a matrice.
> >
> >In the following exampe, if you let RowReduce working as it will, you
> >also reduce the last column
> >(since formally the polynoms in yi are not ZERO) and you cannot get the
> >equations of image of the matrice.
> >
> >In[23]:=n = 5;
> > m = Table[(i + j - 1)^2, {i, n}, {j, n}];
> > v = Table[yi, {i, n}];
> > << "LinearAlgebra`MatrixManipulation`";
> > mv = AppendRows[m, Transpose[{v}]]
> >Out[27]=...
> >
> >In[30]:=RowReduce[mv]
> >
> >Out[31]=...
> >
> >Third and (last) question : >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
> >
> >So I qickly wrote some lines to get my aim.
> >But I am just beginning to program with Mathematica
> >and it is written as with Pascal and other classical languages.
> >Is ther a mean to write such a function using more pattern-matching ?
> >
> >gaussElim[m_, jmax_] :=
> > Module[{},
> > mm = m;
> > nl = Dimensions[mm][[1]];
> > nc = Dimensions[mm][[2]];
> > For[i0 = 1; j0 = 1,
> > i0 <= nl; j0 <= Min[nc, jmax],
> > i0++ ; j0++ ,
> > ip = i0;
> > While[ip < nl && mm[[ip, j0]] == 0, ip++];
> > If[mm[[ip, j0]] == 0, i0--; Continue[]];
> >
> > {mm[[i0]], mm[[ip]]} = {mm[[ip]], mm[[i0]]};
> >
> > For[i = i0 + 1,
> > i <= nl,
> > i++,
> > mm[[i]] = Simplify[mm[[i]] - mm[[i, j0]]/mm[[i0,
> >j0]]*mm[[i0]]]];
> > ];
> > mm]
> >
> >In[29]:= gaussElim[mv, 5]
> >Out[29]= ...
> >
> >
> >Thanks very much for answer,
> >
> >Jack-Michel CORNIL
> >
> >VERSAILLES-FRANCE
> >
--
Murray Eisenberg murray at math.umass.edu
Mathematics & Statistics Dept. phone 413 549-1020 (H)
Univ. of Massachusetts 413 545-2859 (W)
Amherst, MA 01003-4515