netAn4.ma
SingularValue Decomposition for Neural Network Data
Notebook Version: 2.0
Steve Hunka, University of Alberta
steve.hunka@ualberta.ca
Mod. Date: Oct 30/97

General Description of the Use of Singular Value Decomposition

In singular value decomposition of the matrix M into ULW (i.e., M=ULW),  U is a set of left eigenvectors, L a diagonal matrix of eigenvalues, and W the right set of eigenvectors.  A=U(L^0.5) and B=(L^0.5)W (referred to as principal components) are calculated and plotted. When M is the hidden cell responses of the last epoch, A  may provide some understanding of the dimensions along which the hidden cells distinguish the exemplars. An interpretation of A can be made since the rows of these principal components refer to exemplars and their characteristics are known.  The decomposition may be used on hidden cell responses and hidden cell weights. 2D and 3D plots are available. The error of approximating M by successive eigenvalues (roots) and their associated eigenvalues is presented in a 2D plot.
 The matrix A may be rotated orthogonally using the Varimax criterion providing AT where T is the matrix of rotation calculated by Varimax. If AT is formed for the left components, then Inverse[T]B must be used for the corresponding right components so that M=AT (T^-1)B is preserved because TT^-1=Identity matrix. No new information is provided by this rotation, but it may provide for easier interpretation of the left components since axes will tend to go through clusters of exemplars. Examples use data from TCLast written by function bpnn in Neural.nb.

Functions in this Notebook

decompSV   singular value decomposition
refmcol       reflects cols when SS(-)>SS(+)
ntrplot       entry to 3D plotting of points dspts3dv3
dspts3dv3  does the 3D plotting
pairPlot     plots x-y points; Requires <<Grapahics`Graphics`
variMax     Varimax rotation
                 aTan4
                 cos4sin
                 rordcol
                 srtindx
                 (refmcol)
                 
 The following functions are  provided by the Package netUtil.nb and documented in netUdoc.nb:
 pmat2:          to display a matrix 
 genrch:         to generate row/column headings for pmat2 
 rndoff:          round off
 kbin3:           input of numeric values from keyboard
 kbinstrg:       input of string from keyboard
 plotLegend:   plots a color legend of colors used in plotnLists
 plotnLists:    2D line plot

Required Packages

<< "Graphics`Graphics`"<< "netUtil`"

dcompSV[d] does a singular value decomposition of the matrix d (D=ULW). Left and right vectors are rescaled by Sqrt[eigenvalues] and referred to as principal components (PC). The left PC are reflected when SS(neg)>SS(pos) as well as corresponding rows of the right PC. The matrix d may be hidden cell weights (w) or responses (r), and the user is queried for "w" or "r" in order to set the correct row and column titles of output and x-axis label in the plots.  The proportion of sqd error of approximating d by 1,(1,2),(1,2,3),...(1,...n) roots is plotted. The proportion of sqd error in approximating d by the successive roots and their associated eigenvectors is based on  DD'=UL^2U' and Trace(Diag(DD'))=SS(D)=Sum[L^2]. 
2D line-plots are made of the right and left PC for the largest roots selected by the user.
Caution: When d is a matrix of hidden cell weights of order (#hidden cells,#input cells) the right PC will have as many columns as input cells. In the case of d being hidden responses of order (#exemplars,#hidden cells) the left PC will have as many rows as exemplars.
Calls: refmcol, pmat2, genrch, rndoff,kbin3,kbinstrg,plotLegend,plotnLists
Returns: {left pc,right pc} for use in other analyses
Mod. Date: May 9/97

dcompSV[d_] := 
      Module[{svm, lpc, ord, rpc, eig, er, rch, ans, ntp, lxpt, rxpt, lpc1, 
          rpc1}, Print["Singular Value Decomposition"]; 
        Print["Dimensions of input data=", Dimensions[d]]; 
        ans = kbinstrg[
            "Are you using Hidden wts or responses? w/r", {{"w"}, {
                "r"}}]; 
        If["r" == ans, lxpt = "Exemplar"; rxpt = "HC"; 
          Print["Input Data are hidden cell Responses"];, 
          lxpt = "HC"; rxpt = "InC"; 
          Print["Input data are hidden cell Weights"]]; 
        svm = SingularValues[d]; 
        eig = svm[[2]]; {lpc, ord} = 
          refmcol[Transpose[svm[[1]]] . 
              DiagonalMatrix[(Sqrt[eig])]]; 
        Print["\nLeft Principal Components"]; 
        If[ans == "r", rch = genrch["Ex", "PC", Dimensions[lpc]], 
          rch = genrch["HC", "PC", Dimensions[lpc]]]; 
        pmat2[rndoff[lpc, 4], 5, rch]; 
        Print["\nReflection Vector=", ord]; 
        Print["\nEigenvalues (Roots)"]; 
        Print[TableForm[eig, TableSpacing -> {0}, 
            TableHeadings -> Automatic]]; 
        rpc = DiagonalMatrix[(Sqrt[eig])] . 
            svm[[3]]; 
        Print["\nRight Principal Components"]; rpc = rpc ord; 
        If[ans == "r", 
          rch = genrch["PC", "HC", Dimensions[rpc]]; , 
          rch = genrch["PC", "IC", Dimensions[rpc]]]; 
        pmat2[rndoff[rpc, 4], 5, rch]; 
        er = Drop[FoldList[Plus, 0, (eig^2)], 1]; er = 1.  - (er/Last)[er]; 
        Print["\nProportion Sqd-error by Root"]; 
        Print[TableForm[er, TableSpacing -> {0}, 
            TableHeadings -> Automatic]]; 
        ListPlot[er, PlotJoined -> True, 
          AxesLabel -> {"Root", "P(e^2)"}, 
          PlotLabel -> "Prop. Sqd-Error of Approx"]; 
        ans = kbinstrg["Line plot PC? y/n", {{"y"}, {"n"}}]; 
        If["y" == ans, 
          ntp = kbin3["How many of largest PC to plot?", 1]; 
          lpc1 = Transpose[Take[Transpose[lpc], ntp]]; plotLegend[ntp]; 
          plotnLists[lpc1, " Left PC", lxpt, 0, 1, 1]]; 
        rpc1 = Transpose[Take[rpc, ntp]]; 
        plotnLists[rpc1, " Right PC", rxpt, 0, 1, 1]; 
        Return[{lpc, rpc}]; ];

refmcol[m] reflect cols of matrix m when SS neg values>SS pos values; returns reflected cols and vector of 1 and -1 where -1 indicates col reflected.
Called by: variMax, dcompSV
Mod. Date:

refmcol[m_] := 
      Module[{c, r, ssp, ssn, tm, rec}, tm = Transpose[m]; 
        r = Dimensions[tm][[1]]; 
        c = Dimensions[tm][[2]]; 
        rec = Table[1, {r}]; 
        Do[ssp = 0; ssn = 0; 
          Do[If[0 ≤ 
                tm[[j
                    ]][[k
                  ]], 
              ssp += tm[[j
                      ]][[k
                    ]]^2, 
              ssn += tm[[j
                      ]][[k
                    ]]^2], {k, 1, c}]; 
          If[ssn > ssp, 
              tm[[j]] = 
                -1.0 tm[[j]]; 
              rec[[j]] = -1.0];, {
            j, 1, r}]; tm = Transpose[tm]; Return[{tm, rec}]; ];

Test Data for Singular Value Decomposition

We will use the hidden response of the last epoch stored in file TCLast

rprtFdata

\(tclhr = rcvrExpr["TCL953hdr", 19971022214210]; 

    Dimensions[tclhr]\)

Note to SVD which follows

A very good approximation to the matrix of hidden response at the last epoch is obtained by the use of the first 3 eigenvalues and associated eigenvectors.  This can be seen from the plot of the proportion of sqd error of approximation.

{lfthr, rghthr} = dcompSV[tclhr]

Now do the SVD using hidden cell weights

\(tclhw = rcvrExpr["TCL953hwt", 19971022214210]; 

    Dimensions[tclhw]\)

{lft, rght} = dcompSV[tclhw]

2D and 3D Plotting of PC

General Procedure: The 3D plotting of the left principal component (PC) as returned by function dcompSV is made using the function ntrplot[dp]. The user may first rotate the left PC using Varimax and plot the rotated matrix, however the corresponding right PC matrix must be premultiplied by Inverse[T] where T is the varimax transformation (rotation) matrix. The function pairPlotf is used to plot ALL combinations of 2 cols of the matrix in a 2D plot.

ntrplot[dp] is a function to allow entry into 2d and 3d plotting functions; dp is the matrix to be plotted. If there are only 2 cols in dp, a 2D plot is automatically called.  The following options are provided for 3D plots:
a) specification of the combination of 3 cols to plot
b) the size of the point representing the coordinate
c) whether the points are to be numbered
d) whether the mid-plane is opaque or transparent
e) a title for the plot
f) repeat sequence (a) to (d) 
If the user has 3 or more cols and wants 2D plots, call pairPlotf directly and all pairwise plots are made. 
Calls: pairPlotf, dispts3dv3, kbin3,kbinstrg
Mod. Date: May 9/97

ntrplot[dp_] := 
      Module[{nf, rpq, fnumb, ps, tn, tb, ttl}, 
        nf = Dimensions[dp][[2]]; 
        If[2 == nf, pairPlotf[dp]; Return[]]; rpq = "y"; 
        While[rpq == "y", 
          fnumb = kbin3["Type column numbers to plot, eg {1,2,3}", {3}]; 
          If[0 > Min[fnumb] || nf < Max[fnumb], 
            Print["Column number error"]; Return[]]; 
          ps = kbin3["Type point size,eg .02",  .02]; 
          tn = kbinstrg[
              "Points to be numbered? y/n", {{"y"}, {"n"}}]; 
          If[tn == "y", tn = 1, tn = 0]; 
          tb = kbinstrg["Opaque mid-plane? y/n", {{"y"}, {"n"}}]; 
          If[tb == "y", tb = 0, tb = 1]; 
          ttl = InputString["Type short title for plot"]; 
          dispts3dv3[dp, fnumb, ps, tn, tb, ttl]; 
          rpq = kbinstrg[
              "More plots? y/n", {{"y"}, {"n"}}]; ]; ];

dispts3dv3[dp,fn,ps,tn,tb,tl] displays numbered or unnumbered points in a 3D space 
with vertical line from mid-plane to each point. 
dp: points to plot (by rows) matrix
fn:  vector of 3 column indices to plot
ps:  point size;use .01(small) to .05(large)
tn:  1 to display point numbers; else no display
tb:  0 opaque plane, 1 transparent plane
tl:   title for plot
When an opaque mid-plane is used and a point number appears with no associated line, the point is hidden below the mid-plane. To change font add the option DefaultFont->{"font",size}, or {"font-form",size} where form is for example, bold, to the Show operation. A ViewPoint statement is included which the user may wish to change.
Called by: ntrplot

Mod. Date: Feb 12/97

dispts3dv3[dp_, fn_, ps_, tn_, tb_, tl_] := 
      Module[{nf, nt, d, e, f, g, fctaxes, ax, pobj, slctdp, mta, ctr, lo, 
          hi}, nt = 
          Dimensions[dp][[1]]; 
        slctdp = 
          Transpose[
            Transpose[dp][[fn]]]; 
        d = Point/@slctdp; 
        e = Transpose[
            Append[Drop[Transpose[slctdp], -1], 
              Table[0, {Length[slctdp]}]]]; 
        g = Line/@
            Partition[
              Partition[
                Flatten[Transpose[Join[Transpose[slctdp], Transpose[e]]]], 
                3], 2]; lo = Min/@Transpose[slctdp]; 
        hi = Max/@Transpose[slctdp]; 
        f = {{lo[[1]], 
              lo[[2]], 0}, {
              hi[[1]], 
              lo[[2]], 0}, {
              hi[[1]], 
              hi[[2]], 0}, {
              lo[[1]], 
              hi[[2]], 0}}; 
        fctaxes = {{{hi[[1]], 0, 0}, {
                lo[[1]], 0, 0}}, {{0, 
                hi[[2]], 0}, {0, 
                lo[[2]], 0}}, {{0, 0, 
                hi[[3]]}, {0, 0, 
                lo[[3]]}}}; 
        ax = Line/@fctaxes; 
        Which[tb == 0, 
          pobj = Graphics3D[{PointSize[ps], d, g, ax, Polygon[f]}]; , 
          tb == 1, 
          f = {{lo[[1]], 
                lo[[2]], 0}, {
                hi[[1]], 
                lo[[2]], 0}, {
                hi[[1]], 
                hi[[2]], 0}, {
                lo[[1]], 
                hi[[2]], 0}, {
                lo[[1]], 
                lo[[2]], 0}}; 
          pobj = Graphics3D[{PointSize[ps], d, g, ax, Line[f]}]]; mta = {}; 
        If[tn == 1, mta = Text/@slctdp; 
          Do[If[slctdp[[j
                    ]][[3
                  ]] ≥ 0, ctr = -1, 
              ctr = 1]; 
            mta = mta /. 
                  Text[slctdp[[j]]] 
                    -> 
                    Text[j, 
                      slctdp[[j]], {0, 
                        ctr}];, {j, 1, nt}]]; 
        pobj = Join[pobj[[1]], mta]; 
        Show[Graphics3D[pobj], Axes -> True, PlotRange -> All, 
          PlotLabel -> tl, AxesLabel -> fn]; ];

ntrplot[lfthr]

Notes to 3D plot shown above.The character T in its 4 orientations are shown by points 1,2,3, and 4 which appear well separated from the orientations of C and L. (Points 1 and 3, have no vertical lines to the opaque plane because they fall below the surface.) The character C are defined by points 5,6,7 and 8 are also nicely clustered in a region of the 3 space. The character L defined by points 9,10,11, and 12 suggest that it is being treated more like the character C than T, however the first 2 orientations (9 and 10) are clustered in different parts of the space than the last 2 orientations (11,12). For further clarification, see the first 2D plot below which shows that the exemplars of T are distributed along one line, but C and L appear along a different line with the first 2 and last 2 orientations of L appearing at the ends of the continuum for C. Also see the 3D plot after the Varimax rotation and corresponding pairwise plots below.

pairPlotf[p] plots all pairwise columns of the matrix p.
Caution: The number of ways of taking 2 cols out of n cols can be quite large.  
Mod. Date: Apr 7/97

pairPlotf[p_] := 
      Module[{dp, dt, nf, nr}, 
        nf = Dimensions[p][[2]]; 
        nr = Dimensions[p][[1]]; 
        dt = Join[Transpose[p], {Range[nr]}]; 
        Do[Do[dp = 
              Transpose[
                dt[[{j, k, nf + 1}]]]; 
            LabeledListPlot[dp, AxesLabel -> {j, k}];, {k, j + 1, 
              nf}], {j, 1, nf - 1}]; ];

pairPlotf[Transpose[Take[Transpose[lfthr], 3]]]

Varimax Rotation Function

The Varimax rotation is normally used in factor analysis. (See for example, Mulaik, S.A., The Foundations of Factor Analysis, McGraw-Hill, NY, p258, 1972). Varimax iteratively does pairwise orthogonal rotations of the columns of the matrix of coefficients to maximize the variance of the squared coefficients by column. The effect is to place new axes through clusters of points.(This code was developed from a translation of APL code written many years ago, and thus, is probably not very efficient.)

variMax[a] does rotation on matrix a; 
Calls: aTan4, cos4sin, to get sine and cos for transformation matrix
srtindex, reordcol, refmcol, pmat2, genrch, rndoff
Returns: {rotated matrix, transformation matrix}
Mod. Date: May 10/97

variMax[a_] := 
      Module[{c, r, t, h2, a1, at, zr, f, u, v, nm, dn, sc, ss, t1, index, 
          rch, rfindx}, Print["Varimax Rotation"]; 
        c = Dimensions[a][[2]]; 
        r = Dimensions[a][[1]]; 
        t = IdentityMatrix[c]; h2 = Plus@@Transpose[a]^2; 
        Print["\nStarting Communalities"]; 
        Print[TableForm[h2, TableSpacing -> {0}, 
            TableHeadings -> Automatic]]; a1 = (a/(h2^0.5)); 
        at = Transpose[a1]; zr = (1/2) ((c^2) - c); 
        While[zr > 0, 
          Do[Do[v = 
                2 at[[j]] 
                  at[[k]]; 
              u = at[[j]]^2 - 
                  at[[k]]^2; 
              nm = 2 Plus@@(u v) - ((2 Plus@@u Plus@@v)/r); 
              dn = Plus@@((u^2) - (v^2)) - ((Plus@@(u^2) - Plus@@(v^2))/r); 
              sc = aTan4[nm, dn]; 
              If[0 ≠ 
                    sc[[1]], 
                  t1 = IdentityMatrix[c]; 
                  t1[[j, j]] = 
                    sc[[2]]; 
                  t1[[k, k]] = 
                    sc[[2]]; 
                  t1[[k, j]] = 
                    sc[[3]]; 
                  t1[[j, k]] = 
                    -sc[[3]]; 
                  t = t . t1; at = Transpose[Transpose[at] . t1];, 
                  zr = zr - 1];, {k, j + 1, c}]; , {j, 1, c - 1}]; 
          If[zr ≠ 0, zr = (1/2) ((c^2) - c)]; ]; 
        at = Transpose[at] (h2^0.5); index = srtindx[at]; 
        at = reordcol[at, index]; {at, rfindx} = refmcol[at]; 
        Print["\nRotated Varimax Factors"]; 
        rch = genrch["Var", "F", Dimensions[at]]; 
        pmat2[rndoff[at, 4], 5, rch]; ss = Plus@@(at^2); 
        Print["\nFactor SS"]; 
        Print[TableForm[ss, TableSpacing -> {0}, 
            TableHeadings -> Automatic]]; 
        ss = Plus@@Transpose[at]^2; 
        Print["\nCommunalities after Rotation"]; 
        Print[TableForm[ss, TableSpacing -> {0}, 
            TableHeadings -> Automatic]]; 
        Print["\nTransformation Matrix"]; t = reordcol[t, index]; 
        t = t Table[rfindx, {c}]; 
        rch = genrch["A", "F", Dimensions[t]]; 
        pmat2[rndoff[t, 4], 5, rch]; Return[{at, t}]; ];

aTan4 program; n is numerator and d is denominator; called by variMax
Mod. Date: May 9/97

aTan4[n_, d_] := 
      Module[{tan4, cos4, ctn4, sc, tol, cs}, sc = Table[0, {3}]; tol =  .01; 
        If[Abs[n] == Abs[d] And n == 0, Return[sc]]; 
        Which[Abs[n] < Abs[d], tan4 = Abs[n]/Abs[d]; 
          If[tol < tan4, cos4 = (1/(((1 + (tan4^2)))^0.5)); 
              sc = cos4sin[n, d, cos4]; Return[sc], 
              If[d < 0, sc = {1, (1/(2^0.5)), (1/(2^0.5))}; Return[sc], 
                Return[sc]]];, Abs[n] > Abs[d], ctn4 = Abs[d]/Abs[n]; 
          If[ctn4 < tol, cos4 = 0; sc = cos4sin[n, d, cos4]; Return[sc], 
              cos4 = (ctn4/(((1 + (ctn4^2)))^0.5)); sc = cos4sin[n, d, cos4]; 
              Return[sc]; ];, Abs[n] == Abs[d], cos4 = (1/(2^0.5)); 
          sc = cos4sin[n, d, cos4]; Return[sc]; ]; ];

cos4sin Program; n is numerator, d is denominator, cos4 is cos4 theta; called by aTan4; used in variMax. Mod. Date:

cos4sin[n_, d_, cos4_] := 
      Module[{cs, sc, sn}, 
        cs = ((((1/2) (1 + (((((1 + cos4)/2)))^0.5))))^0.5); 
        sn = (((1 - (cs^2)))^0.5); 
        If[d ≤ 0, cs = (cs/(2^0.5)) + (sn/(2^0.5)); 
          sn = (((1 - (cs^2)))^0.5); sc = {1, cs, sn};, 
          sc = {1, cs, sn}]; 
        If[n ≤ 0, sc = {1, cs, -sn}; Return[sc], 
          Return[sc]]; ];

reordcol[a,indx] reorders cols of matrix a by index given in indx
Mod. Date:

reordcol[a_, indx_] := 
      Transpose[
        Transpose[a][[indx]]];

srtindx[a] returns the sorting index based on the decreasing ss of columns of a. The index returned can then be used for the rearrangement of rows and/or columns of a and other matrices. inx can be used as arg in rordrc[] and rordcol[]. Mod. Date:

srtindx[a_] := 
      Module[{ss, inx, t}, ss = Plus@@(a^2); inx = Range[Length[ss]]; 
        Do[If[ss[[i]] < 
              ss[[j]], 
            t = inx[[i]]; 
            inx[[i]] = 
              inx[[j]]; 
            inx[[j]] = t; ], {i, 1, 
            Length[ss] - 1}, {j, i + 1, Length[ss]}]; Return[inx]; ];

Dimensions[lfthr]

Transpose[Take[Transpose[lfthr], 3]]

vrmlfthr = variMax[Transpose[Take[Transpose[lfthr], 3]]]

ntrplot[vrmlfthr[[1]]]

Notes for the 3D plot shown above. The axes have now been positioned through cluster of points. The orientations of the character T have high values on axis 2, while C and L have low values. C and L distributed themselves along axis 1 suggesting that L is being treated more like C than T. The first 2 orientations of L are treated more like the first 2 orientations of C. L is also distinguished along the 3rd dimension with the first two orientations having high values on 3, and last 2 orientations having small negatiave values. (See pairwise plots below for further clarification.)

pairPlotf[vrmlfthr[[1]]]