Menghitung matriks sirkuler ortogonal

8

Dua baris matriks adalah ortogonal jika produk dalamnya sama dengan nol. Panggil sebuah matriks dengan semua baris ortogonal berpasangan, sebuah matriks ortogonal . Sebuah matriks circulant adalah salah satu di mana setiap vektor baris diputar satu elemen ke relatif kanan ke vektor baris sebelumnya. Kami hanya akan tertarik pada matriks di mana entri adalah salah satu -1atau 1.

Tugas

Tulis kode untuk menghitung sebanyak mungkin perbedaan n/2dengan nmatriks ortogonal, edaran mungkin dalam 2 menit (untuk genap n).

Memasukkan

Kode tidak memiliki input. Itu dapat mencoba nilai rata apa pun yang ndisukainya. Misalnya, kode dapat mencoba semua nyang berlipat 4mulai dari yang terkecil dan juga coba n = 2.

Keluaran

Jumlah matriks sirkuler ortogonal yang Anda temukan. Seharusnya juga ada saklar sederhana dalam kode Anda untuk memungkinkannya mengeluarkan matriks sendiri.

Skor

Jumlah matriks sirkuler yang Anda temukan.

Petunjuk

Matriks orthogonal n/2oleh nsirkulant hanya ada ketika nkelipatan 4atau nkurang dari 4.

Contoh matriks circulant ortogonal adalah:

-1  1 -1 -1
-1 -1  1 -1

Kiat untuk pendekatan naif

Pendekatan yang paling naif adalah hanya mengulangi semua matriks yang mungkin. Ini dapat dipercepat dengan menggunakan dua pengamatan berikut.

  • Untuk menguji ortogonalitas dari matriks sirkuler, kita hanya perlu membandingkan setiap baris dengan yang pertama. Ini diimplementasikan dalam kode sampel.

  • Kita dapat mengulangi kata - kata Lyndon dan kemudian jika kita menemukan matriks ortogonal dikalikan dengan jumlah rotasi yang mungkin. Ini adalah ide yang belum diuji sehingga mungkin bermasalah.

Kode sampel

Ini adalah jawaban python yang sangat sederhana dan naif. Saya menjalankannya menggunakan timeout 120.

import itertools
def check_orthogonal(row):
    for i in xrange(1,int(n/2)):
        if (sum(row[j]*row[(j+i) % n] for j in xrange(n)) != 0):
            return False
    return True

counter = 0
for n in xrange(4,33,4):
    for row in itertools.product([-1,1],repeat = n):
        if check_orthogonal(row):
            counter +=1
            print "Counter is ", counter, ". n = ", n

Tes kebenaran

Sebab n = 4,8,12,16,20,24,28, jumlah matriks berbeda yang harus Anda dapatkan adalah 12,40,144,128,80,192,560, masing-masing.

Tingkat kedahsyatan

Menilai dengan kode sampel, saya dengan ini menyajikan dua tingkat kedahsyatan yang dapat dicoba dicapai oleh jawaban apa pun.

  • Kedahsyatan level perak dicapai dengan mendapatkan skor atau 1156 .

  • Tingkat keangkeran emas adalah untuk mendapatkan yang lebih tinggi dari itu.

Bahasa dan perpustakaan

Anda dapat menggunakan bahasa atau perpustakaan apa pun yang Anda suka (yang tidak dirancang untuk tantangan ini). Namun, untuk tujuan penilaian saya akan menjalankan kode Anda di mesin saya jadi tolong berikan instruksi yang jelas untuk bagaimana menjalankannya di Ubuntu.

Mesin Saya Pengaturan waktu akan dijalankan pada mesin saya. Ini adalah instalasi standar Ubuntu pada Prosesor Delapan-Core AMD 8GB 8GB AMD. Ini juga berarti saya harus dapat menjalankan kode Anda.

Memimpin jawaban

  • 332 oleh flawr di Octave

  • 404 oleh RT dalam Python

  • 744 oleh Contoh Solusi menggunakan pypy

  • 1156 oleh Thomas Kwa menggunakan Java . Kedahsyatan tingkat perak!

  • 1588 oleh Reimer Behrends di OCaml . Kedahsyatan tingkat emas!


sumber
Dapatkah Anda benar-benar menemukan satu untuk setiap nyang merupakan kelipatan empat?
flawr
@ flawr Nah, itu pertanyaan yang bagus! Saya tidak tahu tetapi saya ingin tahu.
Dari apa yang saya lihat sekarang (jika skrip saya benar) mereka tampaknya sangat langka. Sejauh yang saya tahu matriks Hadamard Circulant (matriks persegi dengan entri dalam {-1,1}) diperkirakan hanya ada untuk n = 1 dan 4.
flawr
@ flawr Ya pada kedua hal (saya ingin sesuatu yang jarang membuat tantangan menarik). Tapi saya tidak tahu apa-apa yang diketahui tentang n / 2 oleh n matriks yang beredar.
1
Saya punya solusi menggunakan bitmasking yang menurut saya akan bekerja untuk n = 32.
lirtosiast

Jawaban:

5

OCaml, 1588 (n = 36)

Solusi ini menggunakan pendekatan pola bit biasa untuk mewakili vektor -1s dan 1s. Produk skalar seperti biasa dihitung dengan mengambil xor dari vektor dua bit dan mengurangi n / 2. Vektor bersifat ortogonal jika xornya memiliki n / 2 bit yang tepat.

Kata-kata Lyndon sendiri tidak berguna sebagai representasi yang dinormalisasi untuk ini, karena mereka mengecualikan pola apa pun yang merupakan rotasi itu sendiri. Mereka juga relatif mahal untuk dihitung. Oleh karena itu, kode ini menggunakan bentuk normal yang agak sederhana, yang mensyaratkan bahwa urutan nol terpanjang berturut-turut setelah rotasi (atau salah satunya, jika ada kelipatannya) harus menempati bit paling signifikan. Oleh karena itu, bit yang paling tidak signifikan selalu 1.

Perhatikan juga bahwa setiap vektor kandidat harus memiliki setidaknya n / 4 (dan paling banyak 3n / 4). Karena itu, kami hanya mempertimbangkan vektor dengan n / 4 ... n / 2 bit yang ditetapkan, karena kami dapat menurunkan yang lain melalui komplemen dan rotasi (dalam praktiknya, semua vektor tersebut tampaknya memiliki antara n / 2-2 dan n / 2 + 2 yang , tapi itu sepertinya juga sulit dibuktikan).

Kami membangun formulir normal ini dari bit yang paling signifikan, mengamati kendala bahwa setiap sisa nol yang berjalan (disebut "celah" dalam kode) harus mengikuti persyaratan bentuk normal kami. Secara khusus, selama setidaknya satu bit lagi harus ditempatkan, harus ada ruang untuk celah saat ini dan yang lain selain itu setidaknya sebesar celah saat ini atau celah lain yang diamati sejauh ini.

Kami juga mengamati bahwa daftar hasilnya kecil. Oleh karena itu, kami tidak mencoba untuk menghindari duplikat selama proses penemuan, tetapi hanya mencatat hasil dalam set per-pekerja dan menghitung penyatuan set ini di akhir.

Perlu dicatat bahwa biaya runtime algoritma masih tumbuh secara eksponensial dan pada tingkat yang sebanding dengan versi brute-force; apa ini membeli kita pada dasarnya adalah pengurangan oleh faktor konstan, dan datang pada biaya suatu algoritma yang lebih sulit untuk disejajarkan daripada versi brute-force.

Output untuk n hingga 40:

 4: 12
 8: 40
12: 144
16: 128
20: 80
24: 192
28: 560
32: 0
36: 432
40: 640

Program ini ditulis dalam OCaml, untuk dikompilasi dengan:

ocamlopt -inline 100 -nodynlink -o orthcirc unix.cmxa bigarray.cmxa orthcirc.ml

Jalankan ./orthcirc -helpuntuk melihat opsi apa yang didukung oleh program.

Pada arsitektur yang mendukungnya, -fno-PICmungkin menawarkan sedikit peningkatan kinerja tambahan.

Ini ditulis untuk OCaml 4.02.3, tetapi mungkin juga berfungsi dengan versi yang lebih lama (asalkan tidak terlalu lama).


UPDATE: Versi baru ini menawarkan paralelisasi yang lebih baik. Perhatikan bahwa ia menggunakan p * (n/4 + 1)utas pekerja per instance masalah, dan beberapa dari mereka masih akan berjalan jauh lebih pendek daripada yang lain. Nilai pmust a power of 2. Speedup pada 4-8 core minimal (mungkin sekitar 10%), tetapi itu menskala lebih baik ke sejumlah besar core untuk besar n.

let max_n = ref 40
let min_n = ref 4
let seq_mode = ref false
let show_res = ref false
let fanout = ref 8

let bitcount16 n =
  let b2 n = match n land 3 with 0 -> 0 | 1 | 2 -> 1 | _ -> 2 in
  let b4 n = (b2 n) + (b2 (n lsr 2)) in
  let b8 n = (b4 n) + (b4 (n lsr 4)) in
  (b8 n) + (b8 (n lsr 8))

let bitcount_data =
  let open Bigarray in
  let tmp = Array1.create int8_signed c_layout 65536 in
    for i = 0 to 65535 do
      Array1.set tmp i (bitcount16 i)
    done;
    tmp

let bitcount n =
  let open Bigarray in
  let bc n = Array1.unsafe_get bitcount_data (n land 65535) in
  (bc n) + (bc (n lsr 16)) + (bc (n lsr 32)) + (bc (n lsr 48))

module IntSet = Set.Make (struct
  type t = int
  let compare = Pervasives.compare
end)

let worker_results = ref IntSet.empty

let test_row vec row mask n =
  bitcount ((vec lxor (vec lsr row) lxor (vec lsl (n-row))) land mask) * 2 = n

let record vec len n =
  let m = (1 lsl n) - 1 in
  let rec test_orth_circ ?(row=2) vec m n =
    if 2 * row >= n then true
    else if not (test_row vec row m n) then false
    else test_orth_circ ~row:(row+1) vec m n
  in if test_row vec 1 m n &&
        test_orth_circ vec m n then
  begin
    for i = 0 to n - 1 do
      let v = ((vec lsr i) lor (vec lsl (n - i))) land m in
      worker_results := IntSet.add v !worker_results;
      worker_results := IntSet.add (v lxor m) !worker_results
    done
  end

let show vec n =
  for i = 0 to n / 2 - 1 do
    let vec' = (vec lsr i) lor (vec lsl (n - i)) in
    for j = 0 to n-1 do
      match (vec' lsr (n-j)) land 1 with
      | 0 -> Printf.printf "  1"
      | _ -> Printf.printf " -1"
    done; Printf.printf "\n"
  done; Printf.printf "\n"; flush stdout

let rec build_normalized ~prefix ~plen ~gap ~maxgap ~maxlen ~bits ~fn =
  if bits = 0 then
    fn prefix plen maxlen
  else begin
    let room = maxlen - gap - plen - bits in
    if room >= gap && room >= maxgap then begin
      build_normalized
        ~prefix:(prefix lor (1 lsl (plen + gap)))
        ~plen:(plen + gap + 1)
        ~gap:0
        ~maxgap:(if gap > maxgap then gap else maxgap)
        ~maxlen
        ~bits:(bits - 1)
        ~fn;
      if room > gap + 1 && room > maxgap then
        build_normalized ~prefix ~plen ~gap:(gap + 1) ~maxgap ~maxlen ~bits ~fn
    end
  end

let rec log2 = function
| 0 -> -1
| n -> 1 + (log2 (n lsr 1))

let rec test_gap n pat =
  if n land pat = 0 then true
  else if pat land 1 = 0 then test_gap n (pat lsr 1)
  else false

let rec test_gaps n maxlen len =
  let fill k = (1 lsl k) -1 in
  if len = 0 then []
  else if test_gap n ((fill maxlen) lxor (fill (maxlen-len))) then
    len :: (test_gaps n maxlen (len-1))
  else test_gaps n maxlen (len-1)

let rec longest_gap n len =
  List.fold_left max 0 (test_gaps n len len)

let start_search low lowbits maxlen bits fn =
  let bits = bits - (bitcount low) in
  let plen = log2 low + 1 in
  let gap = lowbits - plen in
  let maxgap = longest_gap low lowbits in
  worker_results := IntSet.empty;
  if bits >= 0 then
    build_normalized ~prefix:low ~plen ~gap ~maxgap ~maxlen ~bits ~fn;
  !worker_results

let spawn f x =
  let open Unix in
  let safe_fork () = try fork() with _ -> -1 in
  let input, output = pipe () in
  let pid = if !seq_mode then -1 else safe_fork() in
  match pid with
  | -1 -> (* seq_mode selected or fork() failed *)
     close input; close output; (fun () -> f x)
  | 0 -> (* child process *)
    close input;
    let to_parent = out_channel_of_descr output in
    Marshal.to_channel to_parent (f x) [];
    close_out to_parent; exit 0
  | pid -> (* parent process *)
    close output;
    let from_child = in_channel_of_descr input in
    (fun () -> 
      ignore (waitpid [] pid);
      let result = Marshal.from_channel from_child in
      close_in from_child; result)

let worker1 (n, k) =
  start_search 1 1 n k record

let worker2 (n, k, p) =
  start_search (p * 2 + 1) (log2 !fanout + 1) n k record

let spawn_workers n =
  let queue = Queue.create () in
  if n = 4 || n = 8 then begin
    for i = n / 4 to n / 2 do
      Queue.add (spawn worker1 (n, i)) queue
    done
  end else begin
    for i = n / 2 downto n / 4 do
      for p = 0 to !fanout - 1 do
        Queue.add (spawn worker2 (n, i, p)) queue
      done
    done
  end;
  Queue.fold (fun acc w -> IntSet.union acc (w())) IntSet.empty queue

let main () =
  if !max_n > 60 then begin
    print_endline "error: cannot handle n > 60";
    exit 1
  end;
  min_n := max !min_n 4;
  if bitcount !fanout <> 1 then begin
    print_endline "error: number of threads must be a power of 2";
    exit 1;
  end;
  for n = !min_n to !max_n do
    if n mod 4 = 0 then
    let result = spawn_workers n in
    Printf.printf "%2d: %d\n" n (IntSet.cardinal result);
    if !show_res then
      IntSet.iter (fun v -> show v n) result;
    flush stdout
  done

let () =
  let args =[("-m", Arg.Set_int min_n, "min size of the n by n/2 matrix");
             ("-n", Arg.Set_int max_n, "max size of the n by n/2 matrix");
             ("-p", Arg.Set_int fanout, "parallel fanout");
             ("-seq", Arg.Set seq_mode, "run in single-threaded mode");
             ("-show", Arg.Set show_res, "display list of results") ] in
  let usage = ("Usage: " ^
               (Filename.basename Sys.argv.(0)) ^
               " [-n size] [-seq] [-show]") in
  let error _ = Arg.usage args usage; exit 1 in
  Arg.parse args error usage;
  main ()
Reimer Behrends
sumber
Ini bagus dan selamat datang kembali! Setelah mengatakan itu ... bagaimana dengan jawaban Nim juga? :)
Kode Anda hanya sampai 36: 432 pada mesin saya dalam waktu.
Hah. Ini berjalan ke 40 hanya dalam waktu kurang dari dua menit pada laptop saya dengan prosesor quadcore (2.5GHz Intel Core i7), jadi saya cukup yakin itu akan membuatnya menjadi 40 pada Anda juga. Saya akan memperbarui sesuai. Tentang pertanyaan Anda yang lain, sementara saya memiliki implementasi Nim brute-force, yang satu masih berjalan dua kali lebih lambat (karena bagian brute-force) dan tidak terlalu berbeda dari kode Thomas Kwa (hanya sedikit lebih banyak pengurangan pencarian ruang), jadi Anda tidak kehilangan sesuatu yang mengasyikkan.
Reimer Behrends
Apakah saya benar bahwa kode Anda memerlukan sistem 64 bit? Saya baru saja mengujinya pada bit 32 di mana selalu output 0.
1
Benar, karena menyimpan vektor sebagai int. Saya juga bisa memaksa representasi 64-bit (atau bahkan bilangan bulat besar), tetapi itu akan sangat lambat pada sistem 32-bit.
Reimer Behrends
3

Jawa, 1156 matriks

Ini menggunakan bitmasking yang cukup naif, dan membutuhkan waktu di bawah 15 detik untuk n = 28 pada mesin saya.

Matriks Circulant ditentukan oleh baris pertama mereka. Oleh karena itu, saya mewakili baris pertama dari matriks sebagai vektor bit: 1 dan 0 mewakili -1 dan 1. Dua baris adalah ortogonal ketika jumlah bit yang diset ketika mereka dikoreksi bersama adalah n / 2.

import java.util.Arrays;

class Main {

    static void dispmatrix(long y,int N)
    {
        int[][] arr = new int[N/2][N];
        for(int row=0; row < N/2; row++)
        {
            for(int col=0; col < N; col++)
            {
                arr[row][col] = (int) ((y >>> (N+row-col)) & 1L);
            }
        }
        System.out.println(Arrays.deepToString(arr));
    }

  public static void main(String[] args) {

    int count = 0;
    boolean good;
    boolean display = false;
    long y;
    for(int N=4; N <= 28 ;N += 4)
    {
    long mask = (1L << N) - 1;
    for(long x=0; x < (1L<<N); x++)
    {
        good = true;
        y = x + (x << N);

        for(int row = 1; row < N/2; row++)
        {
            good &= N/2 == Long.bitCount((y ^ (y >>> row)) & mask);
        }

        if(good)
        {
            if(display) dispmatrix(y,N);
            count++;
        }

    }
    System.out.println(count);
    }
  }
}

Saya tidak bisa membuat Eclipse bekerja sekarang, jadi ini diuji pada repl.it.

Berikut adalah jumlah baris pertama yang ortogonal ke baris pertama setelahnya untuk n = 28:

[268435456, 80233200, 23557248, 7060320, 2083424, 640304, 177408, 53088, 14896, 4144, 2128, 1008, 1008, 560]

Optimasi:

  • Karena ortogonalitas tidak berubah dengan rotasi siklik dari kedua baris, kita hanya perlu memeriksa bahwa baris pertama adalah ortogonal bagi yang lain.
  • Daripada secara manual melakukan pergantian siklik N-bit N / 2 kali, saya menyimpan bit yang akan digeser di atas long, dan kemudian menggunakan satu anddengan bit N yang lebih rendah untuk mengekstrak yang diperlukan.

Kemungkinan optimasi lebih lanjut:

  • Hasilkan kata-kata Lyndon. Menurut perhitungan saya ini hanya masuk akal jika kata-kata Lyndon dapat dihasilkan dalam kurang dari ~ 1000 siklus masing-masing.
  • Jika kata-kata Lyndon terlalu lama, kita masih bisa memeriksa vektor bit yang dimulai 00, dan menggunakan rotasi mereka (dan rotasi NOT) ketika kita menemukan matriks ortogonal. Kita dapat melakukan ini karena 0101...01dan 1010...10tidak mungkin baris pertama, dan semua yang lainnya mengandung a 00atau a 11.
  • Cabang mungkin setengah jalan ketika matriks mungkin ortogonal. Saya tidak tahu berapa banyak biaya percabangan, jadi saya harus menguji.
  • Jika cara di atas berhasil, mulailah dengan baris selain dari baris 1?
  • Mungkin ada beberapa cara matematika untuk menghilangkan beberapa kemungkinan. Saya tidak tahu apa yang akan terjadi.
  • Menulis di C ++, tentu saja.
lirtosiast
sumber
Ini bagus. Terima kasih! Saya menantikan beberapa optimasi Anda sehingga kami dapat melihat angka n=36juga.
Oh dan Anda mencapai kedahsyatan tingkat Silver! :)
2

Python (404 matriks pada i5-5300U)

Terutama memposting ini sebagai titik awal untuk meningkatkan orang lain, ini dapat dibersihkan banyak, diparalelkan, dll.

import numpy
import itertools
import time

def findCirculantOrthogonalRows(n, t1, timeLimit):
  validRows = []
  testMatrix = numpy.zeros((n//2, n))
  identityMatrixScaled = n*numpy.identity(n//2)
  outOfTime = False
  for startingRowTuple in itertools.product([1, -1], repeat=n):
     for offset in range(n//2):
       for index in range(n):
         testMatrix[offset][index] = startingRowTuple[(index-offset) % n]

     if(numpy.array_equal(identityMatrixScaled, testMatrix.dot(testMatrix.transpose()))):
      validRows.append(startingRowTuple)

    if(time.clock() - t1 >= timeLimit):
      outOfTime = True
      break

  return (validRows, outOfTime)

n = 4
validFirstRows = []
t1 = time.clock()
timeLimit = 120
fullOutput = True

while(True):
  print('calling with', n)
  (moreRows, outOfTime) = findCirculantOrthogonalRows(n, t1, timeLimit)

  if(len(moreRows) > 0):
    validFirstRows.extend(moreRows)

  if(outOfTime == True):
    break

  n += 4

print('Found', len(validFirstRows), 'circulant orthogonal matrices in', timeLimit, 'seconds')

if(fullOutput):
  counter = 1
  for r in validFirstRows:
    n = len(r)
    matrix = numpy.zeros((n//2, n))
    for offset in range(n//2):
      for index in range(n):
        matrix[offset][index] = r[(index-offset) % n]
    print('matrix #', counter, ':\n', matrix)
    counter += 1
    input()
RT
sumber
Saya menambahkan beberapa kode sampel. Perbaikan pertama yang jelas adalah beralih dari kata-kata Lyndon tapi saya tidak yakin bagaimana kode itu.
2

Matlab / Oktaf, 381/328 matriks

Juga hanya pendekatan naif, mencoba setiap kombinasi yang mungkin.

counter = 0;
%ok: 2,4,8
%none: 
tic
for n=[2,4,8,12,16,20];
    m=n/2;
    N=2^n-1;
    for k=1:N

        k/N; %remove ; in order to see progress
        v=(dec2bin(k,n)-'0')*2-1;               %first row

        z=toeplitz([v(1) fliplr(v(m+2:n))], v); %create circulante matrix
        w = z * z.';
        if norm(w - diag(diag(w))) < eps
            counter = counter+1;
            %n    %remove % in order to see the matrices
            %z
        end
        if toc>120;
            break;
        end
    end
end
counter
cacat
sumber
Dalam oktaf kode mencetak sejumlah besar ke layar yang mungkin memperlambatnya. "ans = ...."
Oh benar, saya lupa menambahkan bahwa: Baris paling menjorok ada ndan z, dua ini dapat dikomentari dengan satu %. Dan kemudian Anda dapat menambahkan ;setelah counter = counter+1dan k/N yang akan menekan output.
flawr