SNAP Library 6.0, User Reference  2020-12-09 16:24:20
SNAP, a general purpose, high performance system for analysis and manipulation of large networks
Static Public Member Functions | Static Private Member Functions | List of all members
TSparseSVD Class Reference

#include <linalg.h>

Static Public Member Functions

static void SimpleLanczos (const TMatrix &Matrix, const int &NumEig, TFltV &EigValV, const bool &DoLocalReortoP=false, const bool &SvdMatrixProductP=false)
 
static void Lanczos (const TMatrix &Matrix, int NumEig, int Iters, const TSpSVDReOrtoType &ReOrtoType, TFltV &EigValV, TFltVV &EigVecVV, const bool &SvdMatrixProductP=false)
 
static void Lanczos2 (const TMatrix &Matrix, int MaxNumEig, int MaxSecs, const TSpSVDReOrtoType &ReOrtoType, TFltV &EigValV, TFltVV &EigVecVV, const bool &SvdMatrixProductP=false)
 
static void SimpleLanczosSVD (const TMatrix &Matrix, const int &CalcSV, TFltV &SngValV, const bool &DoLocalReortoP=false)
 
static void LanczosSVD (const TMatrix &Matrix, int NumSV, int Iters, const TSpSVDReOrtoType &ReOrtoType, TFltV &SgnValV, TFltVV &LeftSgnVecVV, TFltVV &RightSgnVecVV)
 
static void OrtoIterSVD (const TMatrix &Matrix, int NumSV, int IterN, TFltV &SgnValV)
 
static void Project (const TIntFltKdV &Vec, const TFltVV &U, TFltV &ProjVec)
 

Static Private Member Functions

static void MultiplyATA (const TMatrix &Matrix, const TFltVV &Vec, int ColId, TFltV &Result)
 
static void MultiplyATA (const TMatrix &Matrix, const TFltV &Vec, TFltV &Result)
 

Detailed Description

Definition at line 407 of file linalg.h.

Member Function Documentation

void TSparseSVD::Lanczos ( const TMatrix &  Matrix,
int  NumEig,
int  Iters,
const TSpSVDReOrtoType &  ReOrtoType,
TFltV &  EigValV,
TFltVV &  EigVecVV,
const bool &  SvdMatrixProductP = false 
)
static

Definition at line 1134 of file linalg.cpp.

1136  {
1137 
1138  if (SvdMatrixProductP) {
1139  // if this fails, use transposed matrix
1140  IAssert(Matrix.GetRows() >= Matrix.GetCols());
1141  } else {
1142  IAssert(Matrix.GetRows() == Matrix.GetCols());
1143  }
1144  IAssertR(NumEig <= Iters, TStr::Fmt("%d <= %d", NumEig, Iters));
1145 
1146  //if (ReOrtoType == ssotFull) printf("Full reortogonalization\n");
1147  int i, N = Matrix.GetCols(), K = 0; // K - current dimension of T
1148  double t = 0.0, eps = 1e-6; // t - 1-norm of T
1149 
1150  //sequence of Ritz's vectors
1151  TFltVV Q(N, Iters);
1152  double tmp = 1/sqrt((double)N);
1153  for (i = 0; i < N; i++) {
1154  Q(i,0) = tmp;
1155  }
1156  //converget Ritz's vectors
1157  TVec<TFltV> ConvgQV(Iters);
1158  TIntV CountConvgV(Iters);
1159  for (i = 0; i < Iters; i++) CountConvgV[i] = 0;
1160  // const int ConvgTreshold = 50;
1161 
1162  //diagonal and subdiagonal of T
1163  TFltV d(Iters+1), e(Iters+1);
1164  //eigenvectors of T
1165  //TFltVV V;
1166  TFltVV V(Iters, Iters);
1167 
1168  // z - current Lanczos's vector
1169  TFltV z(N), bb(Iters), aa(Iters), y(N);
1170  //printf("svd(%d,%d)...\n", NumEig, Iters);
1171 
1172  if (SvdMatrixProductP) {
1173  // A = Matrix'*Matrix
1174  MultiplyATA(Matrix, Q, 0, z);
1175  } else {
1176  // A = Matrix
1177  Matrix.Multiply(Q, 0, z);
1178  }
1179 
1180  for (int j = 0; j < (Iters-1); j++) {
1181  //printf("%d..\r",j+2);
1182 
1183  //calculates (j+1)-th Lanczos's vector
1184  // aa[j] = <Q(:,j), z>
1185  aa[j] = TLinAlg::DotProduct(Q, j, z);
1186  //printf(" %g -- ", aa[j].Val); //HACK
1187 
1188  TLinAlg::AddVec(-aa[j], Q, j, z);
1189  if (j > 0) {
1190  // z := -aa[j] * Q(:,j) + z
1191  TLinAlg::AddVec(-bb[j-1], Q, j-1, z);
1192 
1193  //reortogonalization
1194  if (ReOrtoType == ssotSelective || ReOrtoType == ssotFull) {
1195  for (i = 0; i <= j; i++) {
1196  // if i-tj vector converget, than we have to ortogonalize against it
1197  if ((ReOrtoType == ssotFull) ||
1198  (bb[j-1] * TFlt::Abs(V(K-1, i)) < eps * t)) {
1199 
1200  ConvgQV[i].Reserve(N,N); CountConvgV[i]++;
1201  TFltV& vec = ConvgQV[i];
1202  //vec = Q * V(:,i)
1203  for (int k = 0; k < N; k++) {
1204  vec[k] = 0.0;
1205  for (int l = 0; l < K; l++) {
1206  vec[k] += Q(k,l) * V(l,i);
1207  }
1208  }
1209  TLinAlg::AddVec(-TLinAlg::DotProduct(ConvgQV[i], z), ConvgQV[i], z ,z);
1210  }
1211  }
1212  }
1213  }
1214 
1215  //adds (j+1)-th Lanczos's vector to Q
1216  bb[j] = TLinAlg::Norm(z);
1217  if (!(bb[j] > 1e-10)) {
1218  printf("Rank of matrix is only %d\n", j+2);
1219  printf("Last singular value is %g\n", bb[j].Val);
1220  break;
1221  }
1222 
1223  for (i = 0; i < N; i++) {
1224  Q(i, j+1) = z[i] / bb[j];
1225  }
1226 
1227  //next Lanzcos vector
1228  if (SvdMatrixProductP) {
1229  // A = Matrix'*Matrix
1230  MultiplyATA(Matrix, Q, j+1, z);
1231  } else {
1232  // A = Matrix
1233  Matrix.Multiply(Q, j+1, z);
1234  }
1235 
1236  //calculate T (K x K matrix)
1237  K = j + 2;
1238  // calculate diagonal
1239  for (i = 1; i < K; i++) d[i] = aa[i-1];
1240  d[K] = TLinAlg::DotProduct(Q, K-1, z);
1241  // calculate subdiagonal
1242  e[1] = 0.0;
1243  for (i = 2; i <= K; i++) e[i] = bb[i-2];
1244 
1245  //calculate 1-norm of T
1246  t = TFlt::GetMx(TFlt::Abs(d[1]) + TFlt::Abs(e[2]), TFlt::Abs(e[K]) + TFlt::Abs(d[K]));
1247  for (i = 2; i < K; i++) {
1248  t = TFlt::GetMx(t, TFlt::Abs(e[i]) + TFlt::Abs(d[i]) + TFlt::Abs(e[i+1]));
1249  }
1250 
1251  //set V to identity matrix
1252  //V.Gen(K,K);
1253  for (i = 0; i < K; i++) {
1254  for (int k = 0; k < K; k++) {
1255  V(i,k) = 0.0;
1256  }
1257  V(i,i) = 1.0;
1258  }
1259 
1260  //eigenvectors of T
1261  TNumericalStuff::EigSymmetricTridiag(d, e, K, V);
1262  }//for
1263 
1264  //printf("\n");
1265 
1266  // Finds NumEig largest eigen values
1267  TFltIntKdV sv(K);
1268  for (i = 0; i < K; i++) {
1269  sv[i].Key = TFlt::Abs(d[i+1]);
1270  sv[i].Dat = i;
1271  }
1272  sv.Sort(false);
1273 
1274  TFltV uu(Matrix.GetRows());
1275  const int FinalNumEig = TInt::GetMn(NumEig, K);
1276  EigValV.Reserve(FinalNumEig,0);
1277  EigVecVV.Gen(Matrix.GetCols(), FinalNumEig);
1278  for (i = 0; i < FinalNumEig; i++) {
1279  //printf("s[%d] = %20.15f\r", i, sv[i].Key.Val);
1280  int ii = sv[i].Dat;
1281  double sigma = d[ii+1].Val;
1282  // calculate singular value
1283  EigValV.Add(sigma);
1284  // calculate i-th right singular vector ( V := Q * W )
1285  TLinAlg::Multiply(Q, V, ii, EigVecVV, i);
1286  }
1287  //printf("done \n");
1288 }
IAssert
#define IAssert(Cond)
Definition: bd.h:262
TLinAlg::Norm
static double Norm(const TFltV &x)
Definition: linalg.cpp:324
IAssertR
#define IAssertR(Cond, Reason)
Definition: bd.h:265
TFlt::GetMx
static double GetMx(const double &Flt1, const double &Flt2)
Definition: dt.h:1444
TInt::GetMn
static int GetMn(const int &Int1, const int &Int2)
Definition: dt.h:1183
TMatrix::GetRows
int GetRows() const
Definition: linalg.h:45
TLinAlg::AddVec
static void AddVec(const double &k, const TFltV &x, const TFltV &y, TFltV &z)
Definition: linalg.cpp:232
TMatrix::GetCols
int GetCols() const
Definition: linalg.h:47
TNumericalStuff::EigSymmetricTridiag
static void EigSymmetricTridiag(TFltV &d, TFltV &e, int n, TFltVV &z)
Definition: linalg.cpp:740
TLinAlg::Multiply
static void Multiply(const TFltVV &A, const TFltV &x, TFltV &y)
Definition: linalg.cpp:428
TFlt::Abs
static double Abs(const double &Flt)
Definition: dt.h:1430
TMatrix::Multiply
void Multiply(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.h:24
TStr::Fmt
static TStr Fmt(const char *FmtStr,...)
Definition: dt.cpp:1599
TLinAlg::DotProduct
static double DotProduct(const TFltV &x, const TFltV &y)
Definition: linalg.cpp:165
TVec::Reserve
void Reserve(const TSizeTy &_MxVals)
Reserves enough memory for the vector to store _MxVals elements.
Definition: ds.h:543
TVec::Add
TSizeTy Add()
Adds a new element at the end of the vector, after its current last element.
Definition: ds.h:602
TVVec::Gen
void Gen(const TSizeTy &_XDim, const TSizeTy &_YDim)
Definition: ds.h:2247
TSparseSVD::MultiplyATA
static void MultiplyATA(const TMatrix &Matrix, const TFltVV &Vec, int ColId, TFltV &Result)
Definition: linalg.cpp:1003
void TSparseSVD::Lanczos2 ( const TMatrix &  Matrix,
int  MaxNumEig,
int  MaxSecs,
const TSpSVDReOrtoType &  ReOrtoType,
TFltV &  EigValV,
TFltVV &  EigVecVV,
const bool &  SvdMatrixProductP = false 
)
static

Definition at line 1290 of file linalg.cpp.

1292  {
1293 
1294  if (SvdMatrixProductP) {
1295  // if this fails, use transposed matrix
1296  IAssert(Matrix.GetRows() >= Matrix.GetCols());
1297  } else {
1298  IAssert(Matrix.GetRows() == Matrix.GetCols());
1299  }
1300  //IAssertR(NumEig <= Iters, TStr::Fmt("%d <= %d", NumEig, Iters));
1301 
1302  //if (ReOrtoType == ssotFull) printf("Full reortogonalization\n");
1303  int i, N = Matrix.GetCols(), K = 0; // K - current dimension of T
1304  double t = 0.0, eps = 1e-6; // t - 1-norm of T
1305 
1306  //sequence of Ritz's vectors
1307  TFltVV Q(N, MaxNumEig);
1308  double tmp = 1/sqrt((double)N);
1309  for (i = 0; i < N; i++)
1310  Q(i,0) = tmp;
1311  //converget Ritz's vectors
1312  TVec<TFltV> ConvgQV(MaxNumEig);
1313  TIntV CountConvgV(MaxNumEig);
1314  for (i = 0; i < MaxNumEig; i++) CountConvgV[i] = 0;
1315  // const int ConvgTreshold = 50;
1316 
1317  //diagonal and subdiagonal of T
1318  TFltV d(MaxNumEig+1), e(MaxNumEig+1);
1319  //eigenvectors of T
1320  //TFltVV V;
1321  TFltVV V(MaxNumEig, MaxNumEig);
1322 
1323  // z - current Lanczos's vector
1324  TFltV z(N), bb(MaxNumEig), aa(MaxNumEig), y(N);
1325  //printf("svd(%d,%d)...\n", NumEig, Iters);
1326 
1327  if (SvdMatrixProductP) {
1328  // A = Matrix'*Matrix
1329  MultiplyATA(Matrix, Q, 0, z);
1330  } else {
1331  // A = Matrix
1332  Matrix.Multiply(Q, 0, z);
1333  }
1334  TExeTm ExeTm;
1335  for (int j = 0; j < (MaxNumEig-1); j++) {
1336  printf("%d [%s]..\r",j+2, ExeTm.GetStr());
1337  if (ExeTm.GetSecs() > MaxSecs) { break; }
1338 
1339  //calculates (j+1)-th Lanczos's vector
1340  // aa[j] = <Q(:,j), z>
1341  aa[j] = TLinAlg::DotProduct(Q, j, z);
1342  //printf(" %g -- ", aa[j].Val); //HACK
1343 
1344  TLinAlg::AddVec(-aa[j], Q, j, z);
1345  if (j > 0) {
1346  // z := -aa[j] * Q(:,j) + z
1347  TLinAlg::AddVec(-bb[j-1], Q, j-1, z);
1348 
1349  //reortogonalization
1350  if (ReOrtoType == ssotSelective || ReOrtoType == ssotFull) {
1351  for (i = 0; i <= j; i++) {
1352  // if i-tj vector converget, than we have to ortogonalize against it
1353  if ((ReOrtoType == ssotFull) ||
1354  (bb[j-1] * TFlt::Abs(V(K-1, i)) < eps * t)) {
1355 
1356  ConvgQV[i].Reserve(N,N); CountConvgV[i]++;
1357  TFltV& vec = ConvgQV[i];
1358  //vec = Q * V(:,i)
1359  for (int k = 0; k < N; k++) {
1360  vec[k] = 0.0;
1361  for (int l = 0; l < K; l++)
1362  vec[k] += Q(k,l) * V(l,i);
1363  }
1364  TLinAlg::AddVec(-TLinAlg::DotProduct(ConvgQV[i], z), ConvgQV[i], z ,z);
1365  }
1366  }
1367  }
1368  }
1369 
1370  //adds (j+1)-th Lanczos's vector to Q
1371  bb[j] = TLinAlg::Norm(z);
1372  if (!(bb[j] > 1e-10)) {
1373  printf("Rank of matrix is only %d\n", j+2);
1374  printf("Last singular value is %g\n", bb[j].Val);
1375  break;
1376  }
1377  for (i = 0; i < N; i++)
1378  Q(i, j+1) = z[i] / bb[j];
1379 
1380  //next Lanzcos vector
1381  if (SvdMatrixProductP) {
1382  // A = Matrix'*Matrix
1383  MultiplyATA(Matrix, Q, j+1, z);
1384  } else {
1385  // A = Matrix
1386  Matrix.Multiply(Q, j+1, z);
1387  }
1388 
1389  //calculate T (K x K matrix)
1390  K = j + 2;
1391  // calculate diagonal
1392  for (i = 1; i < K; i++) d[i] = aa[i-1];
1393  d[K] = TLinAlg::DotProduct(Q, K-1, z);
1394  // calculate subdiagonal
1395  e[1] = 0.0;
1396  for (i = 2; i <= K; i++) e[i] = bb[i-2];
1397 
1398  //calculate 1-norm of T
1399  t = TFlt::GetMx(TFlt::Abs(d[1]) + TFlt::Abs(e[2]), TFlt::Abs(e[K]) + TFlt::Abs(d[K]));
1400  for (i = 2; i < K; i++)
1401  t = TFlt::GetMx(t, TFlt::Abs(e[i]) + TFlt::Abs(d[i]) + TFlt::Abs(e[i+1]));
1402 
1403  //set V to identity matrix
1404  //V.Gen(K,K);
1405  for (i = 0; i < K; i++) {
1406  for (int k = 0; k < K; k++)
1407  V(i,k) = 0.0;
1408  V(i,i) = 1.0;
1409  }
1410 
1411  //eigenvectors of T
1412  TNumericalStuff::EigSymmetricTridiag(d, e, K, V);
1413  }//for
1414  printf("... calc %d.", K);
1415  // Finds NumEig largest eigen values
1416  TFltIntKdV sv(K);
1417  for (i = 0; i < K; i++) {
1418  sv[i].Key = TFlt::Abs(d[i+1]);
1419  sv[i].Dat = i;
1420  }
1421  sv.Sort(false);
1422 
1423  TFltV uu(Matrix.GetRows());
1424  const int FinalNumEig = K; //TInt::GetMn(NumEig, K);
1425  EigValV.Reserve(FinalNumEig,0);
1426  EigVecVV.Gen(Matrix.GetCols(), FinalNumEig);
1427  for (i = 0; i < FinalNumEig; i++) {
1428  //printf("s[%d] = %20.15f\r", i, sv[i].Key.Val);
1429  int ii = sv[i].Dat;
1430  double sigma = d[ii+1].Val;
1431  // calculate singular value
1432  EigValV.Add(sigma);
1433  // calculate i-th right singular vector ( V := Q * W )
1434  TLinAlg::Multiply(Q, V, ii, EigVecVV, i);
1435  }
1436  printf(" done\n");
1437 }
IAssert
#define IAssert(Cond)
Definition: bd.h:262
TLinAlg::Norm
static double Norm(const TFltV &x)
Definition: linalg.cpp:324
TExeTm
Definition: tm.h:355
TFlt::GetMx
static double GetMx(const double &Flt1, const double &Flt2)
Definition: dt.h:1444
TMatrix::GetRows
int GetRows() const
Definition: linalg.h:45
TLinAlg::AddVec
static void AddVec(const double &k, const TFltV &x, const TFltV &y, TFltV &z)
Definition: linalg.cpp:232
TMatrix::GetCols
int GetCols() const
Definition: linalg.h:47
TNumericalStuff::EigSymmetricTridiag
static void EigSymmetricTridiag(TFltV &d, TFltV &e, int n, TFltVV &z)
Definition: linalg.cpp:740
TLinAlg::Multiply
static void Multiply(const TFltVV &A, const TFltV &x, TFltV &y)
Definition: linalg.cpp:428
TFlt::Abs
static double Abs(const double &Flt)
Definition: dt.h:1430
TMatrix::Multiply
void Multiply(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.h:24
TExeTm::GetSecs
double GetSecs() const
Definition: tm.h:366
TLinAlg::DotProduct
static double DotProduct(const TFltV &x, const TFltV &y)
Definition: linalg.cpp:165
TVec::Reserve
void Reserve(const TSizeTy &_MxVals)
Reserves enough memory for the vector to store _MxVals elements.
Definition: ds.h:543
TExeTm::GetStr
const char * GetStr() const
Definition: tm.h:368
TVec::Add
TSizeTy Add()
Adds a new element at the end of the vector, after its current last element.
Definition: ds.h:602
TVVec::Gen
void Gen(const TSizeTy &_XDim, const TSizeTy &_YDim)
Definition: ds.h:2247
TSparseSVD::MultiplyATA
static void MultiplyATA(const TMatrix &Matrix, const TFltVV &Vec, int ColId, TFltV &Result)
Definition: linalg.cpp:1003
void TSparseSVD::LanczosSVD ( const TMatrix &  Matrix,
int  NumSV,
int  Iters,
const TSpSVDReOrtoType &  ReOrtoType,
TFltV &  SgnValV,
TFltVV &  LeftSgnVecVV,
TFltVV &  RightSgnVecVV 
)
static

Definition at line 1454 of file linalg.cpp.

1456  {
1457 
1458  // solve eigen problem for Matrix'*Matrix
1459  Lanczos(Matrix, NumSV, Iters, ReOrtoType, SgnValV, RightSgnVecVV, true);
1460  // calculate left singular vectors and sqrt singular values
1461  const int FinalNumSV = SgnValV.Len();
1462  LeftSgnVecVV.Gen(Matrix.GetRows(), FinalNumSV);
1463  TFltV LeftSgnVecV(Matrix.GetRows());
1464  for (int i = 0; i < FinalNumSV; i++) {
1465  if (SgnValV[i].Val < 0.0) { SgnValV[i] = 0.0; }
1466  const double SgnVal = sqrt(SgnValV[i]);
1467  SgnValV[i] = SgnVal;
1468  // calculate i-th left singular vector ( U := A * V * S^(-1) )
1469  Matrix.Multiply(RightSgnVecVV, i, LeftSgnVecV);
1470  for (int j = 0; j < LeftSgnVecV.Len(); j++) {
1471  LeftSgnVecVV(j,i) = LeftSgnVecV[j] / SgnVal; }
1472  }
1473  //printf("done \n");
1474 }
TSparseSVD::Lanczos
static void Lanczos(const TMatrix &Matrix, int NumEig, int Iters, const TSpSVDReOrtoType &ReOrtoType, TFltV &EigValV, TFltVV &EigVecVV, const bool &SvdMatrixProductP=false)
Definition: linalg.cpp:1134
TVec::Len
TSizeTy Len() const
Returns the number of elements in the vector.
Definition: ds.h:575
TMatrix::GetRows
int GetRows() const
Definition: linalg.h:45
TMatrix::Multiply
void Multiply(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.h:24
TVVec::Gen
void Gen(const TSizeTy &_XDim, const TSizeTy &_YDim)
Definition: ds.h:2247
void TSparseSVD::MultiplyATA ( const TMatrix &  Matrix,
const TFltVV &  Vec,
int  ColId,
TFltV &  Result 
)
staticprivate

Definition at line 1003 of file linalg.cpp.

1004  {
1005  TFltV tmp(Matrix.GetRows());
1006  // tmp = A * Vec(:,ColId)
1007  Matrix.Multiply(Vec, ColId, tmp);
1008  // Vec = A' * tmp
1009  Matrix.MultiplyT(tmp, Result);
1010 }
TMatrix::MultiplyT
void MultiplyT(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.h:34
TMatrix::GetRows
int GetRows() const
Definition: linalg.h:45
TMatrix::Multiply
void Multiply(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.h:24
void TSparseSVD::MultiplyATA ( const TMatrix &  Matrix,
const TFltV &  Vec,
TFltV &  Result 
)
staticprivate

Definition at line 1012 of file linalg.cpp.

1013  {
1014  TFltV tmp(Matrix.GetRows());
1015  // tmp = A * Vec
1016  Matrix.Multiply(Vec, tmp);
1017  // Vec = A' * tmp
1018  Matrix.MultiplyT(tmp, Result);
1019 }
TMatrix::MultiplyT
void MultiplyT(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.h:34
TMatrix::GetRows
int GetRows() const
Definition: linalg.h:45
TMatrix::Multiply
void Multiply(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.h:24
void TSparseSVD::OrtoIterSVD ( const TMatrix &  Matrix,
int  NumSV,
int  IterN,
TFltV &  SgnValV 
)
static

Definition at line 1021 of file linalg.cpp.

1022  {
1023 
1024  int i, j, k;
1025  int N = Matrix.GetCols(), M = NumSV;
1026  TFltVV Q(N, M);
1027 
1028  // Q = rand(N,M)
1029  TRnd rnd;
1030  for (i = 0; i < N; i++) {
1031  for (j = 0; j < M; j++)
1032  Q(i,j) = rnd.GetUniDev();
1033  }
1034 
1035  TFltV tmp(N);
1036  for (int IterC = 0; IterC < IterN; IterC++) {
1037  printf("%d..", IterC);
1038  // Gram-Schmidt
1039  TLinAlg::GS(Q);
1040  // Q = A'*A*Q
1041  for (int ColId = 0; ColId < M; ColId++) {
1042  MultiplyATA(Matrix, Q, ColId, tmp);
1043  for (k = 0; k < N; k++) Q(k,ColId) = tmp[k];
1044  }
1045  }
1046 
1047  SgnValV.Reserve(NumSV,0);
1048  for (i = 0; i < NumSV; i++)
1049  SgnValV.Add(sqrt(TLinAlg::Norm(Q,i)));
1050  TLinAlg::GS(Q);
1051 }
TLinAlg::Norm
static double Norm(const TFltV &x)
Definition: linalg.cpp:324
TRnd
Definition: dt.h:11
TLinAlg::GS
static void GS(TVec< TFltV > &Q)
Definition: linalg.cpp:584
TMatrix::GetCols
int GetCols() const
Definition: linalg.h:47
TRnd::GetUniDev
double GetUniDev()
Definition: dt.h:30
TVec::Reserve
void Reserve(const TSizeTy &_MxVals)
Reserves enough memory for the vector to store _MxVals elements.
Definition: ds.h:543
TVec::Add
TSizeTy Add()
Adds a new element at the end of the vector, after its current last element.
Definition: ds.h:602
TSparseSVD::MultiplyATA
static void MultiplyATA(const TMatrix &Matrix, const TFltVV &Vec, int ColId, TFltV &Result)
Definition: linalg.cpp:1003
void TSparseSVD::Project ( const TIntFltKdV &  Vec,
const TFltVV &  U,
TFltV &  ProjVec 
)
static

Definition at line 1476 of file linalg.cpp.

1476  {
1477  const int m = U.GetCols(); // number of columns
1478 
1479  ProjVec.Gen(m, 0);
1480  for (int j = 0; j < m; j++) {
1481  double x = 0.0;
1482  for (int i = 0; i < Vec.Len(); i++)
1483  x += U(Vec[i].Key, j) * Vec[i].Dat;
1484  ProjVec.Add(x);
1485  }
1486 }
TVec::Len
TSizeTy Len() const
Returns the number of elements in the vector.
Definition: ds.h:575
TVec::Gen
void Gen(const TSizeTy &_Vals)
Constructs a vector (an array) of _Vals elements.
Definition: ds.h:523
TVec::Add
TSizeTy Add()
Adds a new element at the end of the vector, after its current last element.
Definition: ds.h:602
TVVec::GetCols
TSizeTy GetCols() const
Definition: ds.h:2253
void TSparseSVD::SimpleLanczos ( const TMatrix &  Matrix,
const int &  NumEig,
TFltV &  EigValV,
const bool &  DoLocalReortoP = false,
const bool &  SvdMatrixProductP = false 
)
static

Definition at line 1053 of file linalg.cpp.

1055  {
1056 
1057  if (SvdMatrixProductP) {
1058  // if this fails, use transposed matrix
1059  IAssert(Matrix.GetRows() >= Matrix.GetCols());
1060  } else {
1061  IAssert(Matrix.GetRows() == Matrix.GetCols());
1062  }
1063 
1064  const int N = Matrix.GetCols(); // size of matrix
1065  TFltV r(N), v0(N), v1(N); // current vector and 2 previous ones
1066  TFltV alpha(NumEig, 0), beta(NumEig, 0); // diagonal and subdiagonal of T
1067 
1068  printf("Calculating %d eigen-values of %d x %d matrix\n", NumEig, N, N);
1069 
1070  // set starting vector
1071  //TRnd Rnd(0);
1072  for (int i = 0; i < N; i++) {
1073  r[i] = 1/sqrt((double)N); // Rnd.GetNrmDev();
1074  v0[i] = v1[i] = 0.0;
1075  }
1076  beta.Add(TLinAlg::Norm(r));
1077 
1078  for (int j = 0; j < NumEig; j++) {
1079  printf("%d\r", j+1);
1080  // v_j -> v_(j-1)
1081  v0 = v1;
1082  // v_j = (1/beta_(j-1)) * r
1083  TLinAlg::MultiplyScalar(1/beta[j], r, v1);
1084  // r = A*v_j
1085  if (SvdMatrixProductP) {
1086  // A = Matrix'*Matrix
1087  MultiplyATA(Matrix, v1, r);
1088  } else {
1089  // A = Matrix
1090  Matrix.Multiply(v1, r);
1091  }
1092  // r = r - beta_(j-1) * v_(j-1)
1093  TLinAlg::AddVec(-beta[j], v0, r, r);
1094  // alpha_j = vj'*r
1095  alpha.Add(TLinAlg::DotProduct(v1, r));
1096  // r = r - v_j * alpha_j
1097  TLinAlg::AddVec(-alpha[j], v1, r, r);
1098  // reortogonalization if neessary
1099  if (DoLocalReortoP) { } //TODO
1100  // beta_j = ||r||_2
1101  beta.Add(TLinAlg::Norm(r));
1102  // compoute approximatie eigenvalues T_j
1103  // test bounds for convergence
1104  }
1105  printf("\n");
1106 
1107  // prepare matrix T
1108  TFltV d(NumEig + 1), e(NumEig + 1);
1109  d[1] = alpha[0]; d[0] = e[0] = e[1] = 0.0;
1110  for (int i = 1; i < NumEig; i++) {
1111  d[i+1] = alpha[i]; e[i+1] = beta[i]; }
1112  // solve eigne problem for tridiagonal matrix with diag d and subdiag e
1113  TFltVV S(NumEig+1,NumEig+1); // eigen-vectors
1114  TLAMisc::FillIdentity(S); // make it identity
1115  TNumericalStuff::EigSymmetricTridiag(d, e, NumEig, S); // solve
1116  //TLAMisc::PrintTFltV(d, "AllEigV");
1117 
1118  // check convergence
1119  TFltKdV AllEigValV(NumEig, 0);
1120  for (int i = 1; i <= NumEig; i++) {
1121  const double ResidualNorm = TFlt::Abs(S(i-1, NumEig-1) * beta.Last());
1122  if (ResidualNorm < 1e-5)
1123  AllEigValV.Add(TFltKd(TFlt::Abs(d[i]), d[i]));
1124  }
1125 
1126  // prepare results
1127  AllEigValV.Sort(false); EigValV.Gen(NumEig, 0);
1128  for (int i = 0; i < AllEigValV.Len(); i++) {
1129  if (i == 0 || (TFlt::Abs(AllEigValV[i].Dat/AllEigValV[i-1].Dat) < 0.9999))
1130  EigValV.Add(AllEigValV[i].Dat);
1131  }
1132 }
IAssert
#define IAssert(Cond)
Definition: bd.h:262
TLinAlg::Norm
static double Norm(const TFltV &x)
Definition: linalg.cpp:324
TFltKd
TKeyDat< TFlt, TFlt > TFltKd
Definition: ds.h:396
TLinAlg::MultiplyScalar
static void MultiplyScalar(const double &k, const TFltV &x, TFltV &y)
Definition: linalg.cpp:414
TMatrix::GetRows
int GetRows() const
Definition: linalg.h:45
TLinAlg::AddVec
static void AddVec(const double &k, const TFltV &x, const TFltV &y, TFltV &z)
Definition: linalg.cpp:232
TMatrix::GetCols
int GetCols() const
Definition: linalg.h:47
TNumericalStuff::EigSymmetricTridiag
static void EigSymmetricTridiag(TFltV &d, TFltV &e, int n, TFltVV &z)
Definition: linalg.cpp:740
TFlt::Abs
static double Abs(const double &Flt)
Definition: dt.h:1430
TMatrix::Multiply
void Multiply(const TFltVV &B, int ColId, TFltV &Result) const
Definition: linalg.h:24
TLAMisc::FillIdentity
static void FillIdentity(TFltVV &M)
Definition: linalg.cpp:1751
TLinAlg::DotProduct
static double DotProduct(const TFltV &x, const TFltV &y)
Definition: linalg.cpp:165
TVec::Gen
void Gen(const TSizeTy &_Vals)
Constructs a vector (an array) of _Vals elements.
Definition: ds.h:523
TVec::Add
TSizeTy Add()
Adds a new element at the end of the vector, after its current last element.
Definition: ds.h:602
TSparseSVD::MultiplyATA
static void MultiplyATA(const TMatrix &Matrix, const TFltVV &Vec, int ColId, TFltV &Result)
Definition: linalg.cpp:1003
void TSparseSVD::SimpleLanczosSVD ( const TMatrix &  Matrix,
const int &  CalcSV,
TFltV &  SngValV,
const bool &  DoLocalReortoP = false 
)
static

Definition at line 1440 of file linalg.cpp.

1441  {
1442 
1443  SimpleLanczos(Matrix, CalcSV, SngValV, DoLocalReorto, true);
1444  for (int SngValN = 0; SngValN < SngValV.Len(); SngValN++) {
1445  //IAssert(SngValV[SngValN] >= 0.0);
1446  if (SngValV[SngValN] < 0.0) {
1447  printf("bad sng val: %d %g\n", SngValN, SngValV[SngValN]());
1448  SngValV[SngValN] = 0;
1449  }
1450  SngValV[SngValN] = sqrt(SngValV[SngValN].Val);
1451  }
1452 }
TVec::Len
TSizeTy Len() const
Returns the number of elements in the vector.
Definition: ds.h:575
TSparseSVD::SimpleLanczos
static void SimpleLanczos(const TMatrix &Matrix, const int &NumEig, TFltV &EigValV, const bool &DoLocalReortoP=false, const bool &SvdMatrixProductP=false)
Definition: linalg.cpp:1053
The documentation for this class was generated from the following files: