> with(LinearAlgebra) > print(); [&x, Add, Adjoint, BackwardSubstitute, BandMatrix, Basis, BezoutMatrix, BidiagonalForm, BilinearForm, CharacteristicMatrix, CharacteristicPolynomial, Column, ColumnDimension, ColumnOperation, ColumnSpace, CompanionMatrix, ConditionNumber, ConstantMatrix, ConstantVector, Copy, CreatePermutation, CrossProduct, DeleteColumn, DeleteRow, Determinant, Diagonal, DiagonalMatrix, Dimension, Dimensions, DotProduct, EigenConditionNumbers, 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, MatrixExponential, MatrixFunction, MatrixInverse, MatrixMatrixMultiply, MatrixNorm, MatrixPower, MatrixScalarMultiply, MatrixVectorMultiply, MinimalPolynomial, Minor, Modular, Multiply, NoUserValue, Norm, Normalize, NullSpace, OuterProductMatrix, Permanent, Pivot, PopovForm, QRDecomposition, RandomMatrix, RandomVector, Rank, RationalCanonicalForm, 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] > A := Matrix(1..2,1..2,[]) > print(); [4 1] [ ] [3 7] > Eigenvectors(A) > print(); [11 1 (1/2)] [ 1 1 ] [-- + - 21 ] [------------- -------------] [2 2 ] [3 1 (1/2) 3 1 (1/2)] [ ], [- + - 21 - - - 21 ] [11 1 (1/2)] [2 2 2 2 ] [-- - - 21 ] [ ] [2 2 ] [ 1 1 ] # # The main point is that the matrix of eigenvectors is not orthogonal. # # We now will find the Singular Value Decomposition. First we find the orthogonal V matrix. # # > B := MatrixMatrixMultiply(Transpose(A), A) > print(); [25 25] [ ] [25 50] > Eigenvectors(B) > print(); [75 25 (1/2)] [ 25 25 ] [-- + -- 5 ] [-------------- --------------] [2 2 ] [25 25 (1/2) 25 25 (1/2)] [ ], [-- + -- 5 -- - -- 5 ] [75 25 (1/2)] [2 2 2 2 ] [-- - -- 5 ] [ ] [2 2 ] [ 1 1 ] > B2 := Matrix(1..2,1..2,[]) > print(); [ 25 25 ] [-------------- --------------] [25 25 (1/2) 25 25 (1/2)] [-- + -- 5 -- - -- 5 ] [2 2 2 2 ] [ ] [ 1 1 ] # # The problem is that this matrix has orthogonal columns but they are not orthonormal. In Maple, Gram-Schmidt normalizes lengths. # > QRDecomposition(B2) > print(); [ (1/2) ] [ 2 2 ] [ ----------------- - --------------------] [ (1/2) (1/2)] [ / (1/2)\ / (1/2)\ ] [ \5 + 5 / \10 - 2 5 / ] [ ], [ (1/2) / (1/2)\ (1/2) ] [2 \1 + 5 / -1 + 5 ] [------------------- -------------------- ] [ (1/2) (1/2) ] [ / (1/2)\ / (1/2)\ ] [2 \5 + 5 / \10 - 2 5 / ] [ (1/2) ] [ (1/2) / (1/2)\ ] [2 \5 + 5 / ] [------------------------ 0 ] [ (1/2) ] [ 1 + 5 ] [ ] [ (1/2)] [ / (1/2)\ ] [ \10 - 2 5 / ] [ 0 --------------------] [ (1/2) ] [ -1 + 5 ] # > V := Matrix(1..2,1..2,[]) > print(); [ (1/2) ] [ 2 2 ] [ ----------------- - --------------------] [ (1/2) (1/2)] [ / (1/2)\ / (1/2)\ ] [ \5 + 5 / \10 - 2 5 / ] [ ] [ (1/2) / (1/2)\ (1/2) ] [2 \1 + 5 / -1 + 5 ] [------------------- -------------------- ] [ (1/2) (1/2) ] [ / (1/2)\ / (1/2)\ ] [2 \5 + 5 / \10 - 2 5 / ] # > MatrixMatrixMultiply(Transpose(V), V) > print(); [[ 2 [[ / (1/2)\ (1/2) [[ 2 \1 + 5 / 2 2 [[---------- + --------------, - --------------------------------------[[ (1/2) / (1/2)\ (1/2) (1/2)[[5 + 5 2 \5 + 5 / / (1/2)\ / (1/2)\ [[ \5 + 5 / \10 - 2 5 / ] [ (1/2) / (1/2)\ / (1/2)\ ] [ 2 \1 + 5 / \-1 + 5 / ] [ + ----------------------------------------], [ (1/2) (1/2)] [ / (1/2)\ / (1/2)\ ] [ 2 \5 + 5 / \10 - 2 5 / ] [ (1/2) 2 2 - -------------------------------------- (1/2) (1/2) / (1/2)\ / (1/2)\ \5 + 5 / \10 - 2 5 / 2] (1/2) / (1/2)\ / (1/2)\ / (1/2)\ ] 2 \1 + 5 / \-1 + 5 / 4 \-1 + 5 / ] + ----------------------------------------, ------------- + --------------] (1/2) (1/2) (1/2) (1/2) ] / (1/2)\ / (1/2)\ 10 - 2 5 10 - 2 5 ] 2 \5 + 5 / \10 - 2 5 / ] ] ] ] ] ] ] ] > simplify(%) > print(); [1 0] [ ] [0 1] # so we see that V is an orthogonal matrix. # # Now we find the U matrix # > C := MatrixMatrixMultiply(A, Transpose(A)) > print(); [17 19] [ ] [19 58] > Eigenvectors(C) > print(); [75 25 (1/2)] [ 19 19 ] [-- + -- 5 ] [-------------- --------------] [2 2 ] [41 25 (1/2) 41 25 (1/2)] [ ], [-- + -- 5 -- - -- 5 ] [75 25 (1/2)] [2 2 2 2 ] [-- - -- 5 ] [ ] [2 2 ] [ 1 1 ] > C2 := Matrix(1..2,1..2,[]) > print(); [ 19 19 ] [-------------- --------------] [41 25 (1/2) 41 25 (1/2)] [-- + -- 5 -- - -- 5 ] [2 2 2 2 ] [ ] [ 1 1 ] > QRDecomposition(C2) > print(); [ (1/2) ] [ 19 2 38 ] [------------------------ - ------------------------] [ (1/2) (1/2)] [ / (1/2)\ / (1/2)\ ] [5 \125 + 41 5 / 5 \250 - 82 5 / ] [ ], [ (1/2) / (1/2)\ (1/2) ] [ 2 \41 + 25 5 / -41 + 25 5 ] [------------------------- ------------------------ ] [ (1/2) (1/2) ] [ / (1/2)\ / (1/2)\ ] [10 \125 + 41 5 / 5 \250 - 82 5 / ] [ (1/2) ] [ (1/2) / (1/2)\ ] [5 2 \125 + 41 5 / ] [------------------------------- 0 ] [ (1/2) ] [ 41 + 25 5 ] [ ] [ (1/2)] [ / (1/2)\ ] [ 5 \250 - 82 5 / ] [ 0 ------------------------] [ (1/2) ] [ -41 + 25 5 ] # > U := Matrix(1..2,1..2,[]) > print(); [ (1/2) ] [ 19 2 38 ] [------------------------ - ------------------------] [ (1/2) (1/2)] [ / (1/2)\ / (1/2)\ ] [5 \125 + 41 5 / 5 \250 - 82 5 / ] [ ] [ (1/2) / (1/2)\ (1/2) ] [ 2 \41 + 25 5 / -41 + 25 5 ] [------------------------- ------------------------ ] [ (1/2) (1/2) ] [ / (1/2)\ / (1/2)\ ] [10 \125 + 41 5 / 5 \250 - 82 5 / ] # # > MatrixMatrixMultiply(Transpose(U), U) > print(); [[ 2 [[ / (1/2)\ [[ 722 \41 + 25 5 / [[-------------------- + --------------------, [[ / (1/2)\ / (1/2)\ [[25 \125 + 41 5 / 50 \125 + 41 5 / [[ (1/2) 722 2 - ------------------------------------------------ (1/2) (1/2) / (1/2)\ / (1/2)\ 25 \125 + 41 5 / \250 - 82 5 / ] [ (1/2) / (1/2)\ / (1/2)\ ] [ 2 \41 + 25 5 / \-41 + 25 5 / ] [ + ------------------------------------------------], [ (1/2) (1/2)] [ / (1/2)\ / (1/2)\ ] [ 50 \125 + 41 5 / \250 - 82 5 / ] [ (1/2) 722 2 - ------------------------------------------------ (1/2) (1/2) / (1/2)\ / (1/2)\ 25 \125 + 41 5 / \250 - 82 5 / (1/2) / (1/2)\ / (1/2)\ 2 \41 + 25 5 / \-41 + 25 5 / + ------------------------------------------------, (1/2) (1/2) / (1/2)\ / (1/2)\ 50 \125 + 41 5 / \250 - 82 5 / 2 ]] / (1/2)\ ]] 1444 \-41 + 25 5 / ]] -------------------- + --------------------]] / (1/2)\ / (1/2)\]] 25 \250 - 82 5 / 25 \250 - 82 5 /]] ]] > simplify(%) > print(); [1 0] [ ] [0 1] # so we see that U is also an orthogonal matrix. # so we see that U is an orthogonal matrix. # # Now we find the diagonal Sigma matrix. # > S := simplify(MatrixMatrixMultiply(MatrixMatrixMultiply(Transpose(U), A), V)) > print(); [5 5 (1/2) ] [- + - 5 0 ] [2 2 ] [ ] [ 5 5 (1/2)] [ 0 - - + - 5 ] [ 2 2 ] > evalf(S) > print(); [8.090169942 0.] [ ] [ 0. 3.090169942] # Now we verify that A = U S V inverse as it is supposed to be. # > simplify(MatrixMatrixMultiply(MatrixMatrixMultiply(U, S), Transpose(V))) > print(); [4 1] [ ] [3 7] # # By the way, singular values are built into Maple > > SingularValues(A) > print(); [ 5 5 (1/2) ] [ - + - 5 ] [ 2 2 ] [ ] [ 5 5 (1/2)] [- - + - 5 ] [ 2 2 ] # # Let's use the power of Maple to calculate some large SVD values. We will start with a random matrix. > A := Vector[Column](1..4,[]) > print(); [ 15 x 15 Matrix ] [ Data Type: anything ] [ Storage: rectangular ] [ Order: Fortran_order ] > evalm(A) > print(); [[28, 33, -48, -49, 10, 85, 89, -13, -6, -68, -81, 80, -57, -46, -15], [33, -65, 28, 52, -99, 27, -12, -96, -49, 66, 97, -89, -31, -60, -89], [-73, -82, 10, -18, -50, 23, 11, -5, -33, 55, -28, 24, 47, -86, 66], [77, 12, 88, -52, -11, 70, -92, -53, -61, 70, -73, 67, 31, 64, 77], [-9, 47, -53, 26, 79, -69, -10, -42, -11, 91, -56, 50, 63, -40, -80], [74, -6, -20, 37, 61, 20, -89, -10, -85, 85, 55, 30, -90, -43, -19], [83, 2, 62, 31, 33, -7, 51, 73, -81, 96, -60, -28, -11, 23, -62], [57, 40, -70, -15, -5, 65, -61, -53, 97, 71, -60, 82, 75, 68, 81], [-97, 20, 55, -46, 69, -86, -60, -55, 10, -20, -76, -68, 33, -2, 22], [-13, -64, -17, 97, 19, 98, 66, 58, -72, 34, -12, 86, 27, 80, 50], [10, 71, 24, -22, -87, 19, 29, 44, 2, 23, 69, 95, 42, -69, 78], [-65, -87, -55, 73, -75, -42, 63, 29, 74, -53, -78, 32, -86, 7, -8], [-80, -96, 65, 90, 50, -5, -48, 10, -83, 65, 81, -41, -77, 86, -90], [-39, -80, 81, 33, -98, -96, -79, -1, -66, -54, 11, 77, -48, -19, -81], [25, -27, -52, -25, 80, 73, 48, 79, 97, 40, 39, 61, -27, -53, -43]] # > SingularValues(A) > print(); [ 1 .. 15 Vector[column] ] [ Data Type: anything ] [ Storage: rectangular ] [ Order: Fortran_order ] > evalm(%) > print(); [9.483260297, 28.47335747, 70.18766891, 109.8491398, 132.0895610, 164.5425943, 185.2450883, 201.4345853, 212.1696226, 220.9141379, 269.0523091, 285.9768444, 326.2782253, 349.2446599, 422.4422120] # Clearly there is not much pattern here and so no dominating singular values, as one might expect; it is hard to compactify noise.. # Try one with more pattern. > A := Matrix(1..10,1..10,[]) > print(); [1 0 1 1 1 1 1 0 0 0] [ ] [0 1 1 1 1 1 1 1 0 0] [ ] [1 0 0 0 0 0 0 0 1 0] [ ] [1 0 3 0 5 0 3 0 1 0] [ ] [1 0 0 0 5 0 0 0 1 0] [ ] [1 0 4 5 5 5 4 0 1 0] [ ] [1 0 0 4 4 4 0 0 1 0] [ ] [0 1 0 0 0 0 0 1 0 0] [ ] [0 0 1 1 1 1 1 0 0 0] [ ] [0 0 1 1 1 1 1 0 0 0] > evalf(SingularValues(A)) > print(); [ 0.] [ ] [ 0.] [ ] [ 0.] [ ] [ 0.] [ ] [0.6552574473] [ ] [ 1.480210675] [ ] [ 1.992485358] [ ] [ 4.353896431] [ ] [ 5.514609483] [ ] [ 14.21415788] # Here the main singular value predominates and the first three singular values prresumably carry most of the pattern. Unfortunately, I do not know how to get Maple to tell me what the associated U and V matrices are, so it is not easy for me to reconstruct the pattern with the first three singular values. . .