> read("/afs/eos.ncsu.edu/users/k/kaltofen/www/courses/LinAlgebra/Maple/
> LsqPkg/InitPkg.mpl");

  libname := /afs/eos.ncsu.edu/users/k/kaltofen/www/courses/LinAlg\

        ebra/Maple, "/afs/bp.ncsu.edu/dist/maple70/lib"

> with(LsqPkg);

                            [Gram_Schmidt]

> A := <<1,1,1,1>|<1,2,3,4>|<1,4,9,16>>;

                              [1    1     1]
                              [            ]
                              [1    2     4]
                         A := [            ]
                              [1    3     9]
                              [            ]
                              [1    4    16]

> with(LinearAlgebra);

  [Add, Adjoint, BackwardSubstitute, BandMatrix, Basis, BezoutMatrix,

        BidiagonalForm, BilinearForm, CharacteristicMatrix,

        CharacteristicPolynomial, Column, ColumnDimension,

        ColumnOperation, ColumnSpace, CompanionMatrix,

        ConditionNumber, ConstantMatrix, ConstantVector,

        CreatePermutation, CrossProduct, DeleteColumn, DeleteRow,

        Determinant, DiagonalMatrix, Dimension, Dimensions,

        DotProduct, Eigenvalues, Eigenvectors, Equal,

        ForwardSubstitute, FrobeniusForm, GaussianElimination,

        GenerateEquations, GenerateMatrix, GetResultDataType,

        GetResultShape, GivensRotationMatrix, GramSchmidt,

        HankelMatrix, HermiteForm, HermitianTranspose, HessenbergForm,

        HilbertMatrix, HouseholderMatrix, IdentityMatrix,

        IntersectionBasis, IsDefinite, IsOrthogonal, IsSimilar,

        IsUnitary, JordanBlockMatrix, JordanForm, LA_Main,

        LUDecomposition, LeastSquares, LinearSolve, Map, Map2,

        MatrixAdd, MatrixInverse, MatrixMatrixMultiply, MatrixNorm,

        MatrixScalarMultiply, MatrixVectorMultiply, MinimalPolynomial,

        Minor, Multiply, NoUserValue, Norm, Normalize, NullSpace,

        OuterProductMatrix, Permanent, Pivot, QRDecomposition,

        RandomMatrix, RandomVector, Rank, ReducedRowEchelonForm, Row,

        RowDimension, RowOperation, RowSpace, ScalarMatrix,

        ScalarMultiply, ScalarVector, SchurForm, SingularValues,

        SmithForm, SubMatrix, SubVector, SumBasis, SylvesterMatrix,

        ToeplitzMatrix, Trace, Transpose, TridiagonalForm, UnitVector,

        VandermondeMatrix, VectorAdd, VectorAngle,

        VectorMatrixMultiply, VectorNorm, VectorScalarMultiply,

        ZeroMatrix, ZeroVector, Zip]

> ?QRDecomposition
> QRDecomposition(A);

     [1/2    - 3/10 sqrt(5)    1/2 ]
     [                             ]  [2       5          15    ]
     [1/2    - 1/10 sqrt(5)    -1/2]  [                         ]
     [                             ], [0    sqrt(5)    5 sqrt(5)]
     [1/2     1/10 sqrt(5)     -1/2]  [                         ]
     [                             ]  [0       0           2    ]
     [1/2     3/10 sqrt(5)     1/2 ]

> ?Gram_Schmidt
> printlevel:=3:
> Q := Gram_Schmidt(A,'R');

                                     [1]
                                     [ ]
                                     [1]
                              v[1]=, [ ]
                                     [1]
                                     [ ]
                                     [1]


                                     [1]
                                     [ ]
                                     [1]
                              u[1]=, [ ]
                                     [1]
                                     [ ]
                                     [1]


                           ||u[1]||^2 = , 4


                                     [1]
                                     [ ]
                                     [2]
                              v[2]=, [ ]
                                     [3]
                                     [ ]
                                     [4]


                                   [-3/2]
                                   [    ]
                                   [-1/2]
                            u[2]=, [    ]
                                   [1/2 ]
                                   [    ]
                                   [3/2 ]


                           <u[1], u[2]>=, 0


                           ||u[2]||^2 = , 5


                                    [ 1]
                                    [  ]
                                    [ 4]
                             v[3]=, [  ]
                                    [ 9]
                                    [  ]
                                    [16]


                                    [ 1]
                                    [  ]
                                    [-1]
                             u[3]=, [  ]
                                    [-1]
                                    [  ]
                                    [ 1]


                           <u[1], u[3]>=, 0


                           <u[2], u[3]>=, 0


                           ||u[3]||^2 = , 4


                             [1    -3/2     1]
                             [               ]
                             [1    -1/2    -1]
                        Q := [               ]
                             [1    1/2     -1]
                             [               ]
                             [1    3/2      1]

> R;

                          [1    5/2    15/2]
                          [                ]
                          [0     1      5  ]
                          [                ]
                          [0     0      1  ]

> Q.R;

                            [1    1     1]
                            [            ]
                            [1    2     4]
                            [            ]
                            [1    3     9]
                            [            ]
                            [1    4    16]

> b := <5,7,29,63>;

                                   [ 5]
                                   [  ]
                                   [ 7]
                              b := [  ]
                                   [29]
                                   [  ]
                                   [63]

> unprotect(D):
> D := Transpose(Q) . Q;

                               [4    0    0]
                               [           ]
                          D := [0    5    0]
                               [           ]
                               [0    0    4]

> D^(-1);

                         [1/4     0      0 ]
                         [                 ]
                         [ 0     1/5     0 ]
                         [                 ]
                         [ 0      0     1/4]

> btilde:=D^(-1) . Transpose(Q) . b;

                                     [ 26 ]
                                     [    ]
                           btilde := [98/5]
                                     [    ]
                                     [ 8  ]

> R;

                          [1    5/2    15/2]
                          [                ]
                          [0     1      5  ]
                          [                ]
                          [0     0      1  ]

> xhat:=R^(-1) . btilde;

                                   [  17  ]
                                   [      ]
                           xhat := [-102/5]
                                   [      ]
                                   [  8   ]

> LinearSolve(Transpose(A) . A, Transpose(A) . b);

                               [  17  ]
                               [      ]
                               [-102/5]
                               [      ]
                               [  8   ]

>