Численный анализ многокубитных систем

[Pipa]

Чаще всего будут использоваться тензорные произведения типа ρ×I и I×ρ, где I – единичная матрица, а ρ – комплексная матрица плотности (I и ρ разной размерности).
Здесь хорошо бы подумать, как ускорить алгоритм и не множить две полных матрицы.

    Готова функция тензорного (прямого) произведения матрицы на единичную справа! Операция умножения в ней отсутствует.
 

Код:

void DirectMulOneLeftR( double **mat1, int n1, int n2, double **mat)
// прямое умножение действительной матрицы на единичную СЛЕВА (M x I)

// **mat1 = 1-я матрица
// n1 = размерность 1-ой матрицы
// n2 = размерность единичной матрицы
// **mat = матрица результата, размерности n1*n2
{
  for( int i=0, ii=0; i < n1; i++, ii+=n2)
  { for( int j=0, jj=0; j < n1; j++, jj+=n2)
    { double x = mat1[i][j];
      for( int k=0; k < n2; k++)
        for( int l=0; l < n2; l++)
          mat[ii+k][jj+l] = (k==l)? x : 0;
    }      
  }
}

Тестирование выглядит так:

Код:

а[][] =
1  2
3  4

DirectMulOneLeftR( a, 2, 3, c);

c[][] =
  1       0        0      2      0       0
  0       1        0      0      2       0
  0       0        1      0      0       2
  3       0        0      4      0       0
  0       3        0      0      4       0
  0       0        3      0      0       4

[С.И. Доронин]

    Хотелось бы уточнить какой (максимальный) размер предполагает у матриц-сомножителей, а также у продукта тензорного произведения (до редукции).

Максимальный размер любых наших матриц не будет превышать размер 8- (или 9) кубитной матрицы. Т.е. той, на которой мы остановимся в качестве максимальной. У сомножителей при тензорном произведении тоже размер будет такой, чтобы результирующая матрица не превышала 8-9 кубитов.

[Pipa]

    Готова функция тензорного (прямого) произведения матрицы на единичную справа! Операция умножения в ней отсутствует.

Код:

void DirectMulOneRightR( double **mat1, int n1, int n2, double **mat)
// прямое умножение действительной матрицы на единичную СПРАВА (I x M)
// **mat1 = 1-я матрица
// n1 = размерность 1-ой матрицы
// n2 = размерность единичной матрицы
// **mat = матрица результата, размерности n1*n2
{
  for( int i=0, ii=0; i < n2; i++, ii+=n1)
    for( int j=0, jj=0; j < n2; j++, jj+=n1)
      for( int k=0; k < n1; k++)
        for( int l=0; l < n1; l++)
            mat[ii+k][jj+l] = (i==j)? mat1[k][l] : 0;
}

Тестирование выглядит так:

Код:

а[][] =
1  2
3  4

DirectMulOneRightR( a, 2, 3, c);

c[][] =
  1       2        0      0      0       0
  3       4        0      0      0       0
  0       0        1      2      0       0
  0       0        3      4      0       0
  0       0        0      0      1       2
  0       0        0      0      3       4

[С.И. Доронин]

Интересно, а какой у них выигрыш по времени по сравнению с произведением двух полных матриц?

[Pipa]

    Одиночные умножения, в отличии от векторных, потребляют очень мало времени. Прямое (тензорное) умножение матриц можно считать мгновенной процедурой (как в случае заполненных матриц, так и с единичной матрицей в качестве сомножителя). Время выполнения таких операций пренебрежимо мало, и находится на уровне погрешностей воспроизведения.

[С.И. Доронин]

Отлично! Хоть одной проблемой меньше , а то я боюсь, что второй блок займет больше времени вычислений, чем первый. Не понятно еще, сколько времени «съест» редукция.

[Pipa]

   На редукцию мне даже страшно смотреть.

[Pipa]

Теперь насчет общего алгоритма редуцирования. Поступаем следующим образом, если мы редуцируем, скажем пятикубитную систему ABCDE и оставляем кубиты AD, проводя редукцию по всем остальным (по ВСЕ), то перебираем все элементы 5-кубитной матрицы плотности и по их индексам в двоичной форме смотрим, что стоит на позициях ВСЕ, т.е. на второй, третьей и пятой позиции в проекторе. Например,
|11111><11111|
Я выделил то, на что нужно смотреть. Если в проекторе слева и справа есть полное совпадение, то этот элемент остается в редуцированной матрице AD и идет, в данном случае к элементам |11><11| (то, что не выделено жирным), таких элементов может получаться несколько, и они будут суммироваться...

    Можно ли понимать эти ваши слова так, что какой бы большой ни была заданная матрица плотности – хоть ABCDEF…XYZ (от A до Z включительно), то редукция ее по AD не превзойдет своим размером 4x4 (т.е. матрицы, изначально построенной на паре AD)? Ведь иначе будет невозможно то "суммирование", про которое вы пишите (т.к. суммирование определено только на матрицах одного размера).
    Вопрос не праздный, поскольку верхний размер редуцированной матрицы желательно знать наперед (из-за того, что практически вырезать столбцы и строки из большой матрицы довольно накладно по времени).

[С.И. Доронин]

Пипа, я не совсем понял вопрос, на выходе, т.е. в результате редуцирования, получается одна матрица, например 4*4, только это не редукция по AD, а наоборот, редукция по ВСЕ, когда остается AD. Исходная матрица, т.е. та, которая редуцируется, не будет превышать наши максимальные 8-9 кубитов.
Варианты редукции могут быть разные, т.е. из одной матрицы, напр. ABCDE (32*32) в результате редукции могут получаться матрицы меньшего размера, напр. А (2*2), AB (4*4), ABC (8*8 ), ABCD (16*16), и аналогичные матрицы для всех других «буквенных» комбинаций.
Т.е. из исходной матрицы получаются все возможные матрицы меньшего размера.

[Pipa]

    Конечно же я неверно выразилась. Под "редукцией по AD" я имела ввиду, что AD остается, а остальное редуцируется.
    Размер исходной матрицы меня не интересует, т.к. ее размер задан при входе в процедуру. Интерес представляет размер продукта. Должна ли сама процедура его вычислять, или же в качестве приемника для редуцированной матрицы может быть подставлена матрица с заранее установленными габаритами?
    Поскольку число "остающихся" от редукции кубитов определяет верхний размер редуцированной матрицы, то это сильно облегчает алгоритм.
    Остается лишь вопрос о необходимости дальнейшего урезания продукта уже по признаку полностью нулевых строк и столбцов. Для каких целей нужна затем редуцированная матрица? Если только для нахождения собственных значений, то урезание нулевых строк/столцов можно не делать - на результат это не повлияет (нулевые соб.значения в расчет не берем).

[С.И. Доронин]

Размер выходной матрицы определяется по входным параметрам. Если в параметрах указаны номера остающихся кубитов (напр. 1 и 4), или их символьное обозначение (AD) то этим и задается размер результирующей матрицы 4*4. Т.е. это то, что говорите Вы:

Поскольку число "остающихся" от редукции кубитов определяет верхний размер релуцированной матрицы, то это сильно облегчает алгоритм.

Для каких целей нужна затем редуцированная матрица? Если только для нахождения собственных значений, то урезание нулевых строк/столцов можно не делать - на результат это не повлияет (нулевые соб.значения в расчет не берем).

Перед тем как находить собственные значения, с этой матрицей нужно будет делать тензорное умножение, поэтому нужен ее настоящий размер.

[Pipa]

Перед тем как находить собственные значения, с этой матрицей нужно будет делать тензорное умножение, поэтому нужен ее настоящий размер.

    Это меняет дело. А то я думала, что сначала делаем тензорное умножение, а лишь потом редуцируем произведение.

[Pipa]

   Впрочем, если делается тензорное умножение на единичную СПРАВА, то собственные значения получатся как объединение собственных значений каждого из сомножителей (поскольку блоки будут выстроены по диагонали). А следовательно, собственные значения исходной матрицы всего лишь дуплицируются столько раз, каков порядок единичной матрицы.
   Что касается умножения на единичную СЛЕВА, то доказать этого не могу, но все числовые примеры дают на MatLab'е тот же результат - набор собственных значений исходной матрицы "размножается" во столько раз, каков порядок единичной матрицы.
   В обоих случаях - как при левом, так и при правом тензорном умножении на единичную матрицу - собственные значения продукта с легкостью угадываются без проведения вычислений в лоб.
   Может быть в этом случае тензорное произведение не надо вычислять?

[С.И. Доронин]

Пипа, перед вычислением собственных значений еще некоторые действия будут делаться. Здесь на форуме я не буду об этом говорить, мы уже вплотную подошли к нашему ноу-хау . Тензорное произведение нужно делать, собственные значения будут находиться для другой матрицы.

[Pipa]

Процедура редуцирования матрицы плотности
(1-ая примерка)

Код:

bool Reduce( double **Re, double **Im, int n1, int offmask,
  double **redRe, double **redIm, int n2)
// input: complex matrix Re:Im, size n1, reduced bites
// output: complex matrix redRe:redIm, size n2
{
// нередуцируемые кубиты
  int onmask = (n1-1)^offmask;  // onmask=9(01001), offmask=22(10110)

Здесь параметр offmask кодирует редуцируемые кубиты. По двоичной системе: A=1, B=2, C=4, D=8, E=16. Редукция системы ABCDE по BCE будет соответствовать коду BCE=2(B)+4(C)+16(E)=22
onmask – комплиментарное к offmask дополнение (кубиты не подвергающиеся редукции). Вычисляется внутри процедуры. Для данного примера onmask = 9 = 1(A)+8(D).

Код:

// массив индексов для переноса
  int *index = new int[n1];
  int n = 0;
  index[0] = 0;
  for( int i=0; i < n1; i++)
    index[i] = ( (i & onmask) && !(i & offmask) )? ++n : 0;
  if( ++n != n2)  // проверка наличия места
  { print( "Bad size: n=%d, n2=%d", n, n2);
    return false;
  }

Массив index[] служит для пересчета номеров столбцов/строк НЕредуцируемых элементов на новые номера в редуцированной матрице. Не все элементы массива используются, но так быстрее будет потом пересчитывать.
Для данного примера:
index[0] = 0; (A-,D-)
index[1] = 1; (A+,D-)
index[8] = 2; (A-,D+)
index[9] = 3; (A+,D+)
остальные элементы  - мусор.

Код:

// обнуляем результат
  for( int i=0; i < n2; i++)
    for( int j=0; j < n2; j++)
      redRe[i][j] = redIm[i][j] = 0;

// добавляем в результат
  for( int i=0; i < n1; i++)
  { int on1 = i & onmask;  // строка: остающиеся
    int off1 = i & offmask;  // строка: редуцируемые
    for( int j=0; j < n1; j++)
    { int on2 = j & onmask;  // столбец: остающиеся
      int off2 = j & offmask;  // столбец: редуцируемые
      if( off1 == off2)  // редуцируемые в столбце и строке совпадают
      { int n = index[on1];  // строка куда
        int m = index[on2];  // столбец куда
        redRe[n][m] += Re[i][j];
        redIm[n][m] += Im[i][j];
      }
    }
  }  

  delete[] index;
  return true;
}

Конец процедуры.

По операциям логического AND с маской offmask выделяются редуцируемые биты (off), а с маской onmask – остающиеся (on). Первые используются для теста "сомнительных элементов" – пропускать их (при различных off) или добавлять (при совпадении off). Вторые используются для расчета места, куда добавлять при совпадении.
Для "базовых" кубитов (AD) исключение не делается – они пропускаются на том же правиле, поскольку их нулевая off должна пропустить их в результат.

Удаление нулевых строк/столбцов результата пока в процедуре отсутствует. Чтобы это сделать необходимо будет завести спец. класс матриц, у которых можно будет выборочно удалять столбцы и строки. Попутно это уменьшит во всей программе число параметров в функциях, поскольку один объект будет хранить в себе действительную, мнимую части и габариты матрицы.

Уверенности в том, что это заработает, у меня нет. Не знаю даже как это проверять. Жду совета.