Schneller Weg zur Berechnung (1: N) '* (1: N)

MATLABs Matrixmultiplikation ist im Allgemeinen ziemlich schnell, aber hier sind ein paar Wege, um nur die obere Dreiecksmatrix zu erhalten. Sie sind langsamer als naiv die v'*v zu berechnen (oder einen MEX-Wrapper zu verwenden, der den passenderen Namen symmetrische Rang k Update-Funktion in BLAS , nicht überraschend!). Wie auch immer, hier sind ein paar MATLAB-only Lösungen:

Die erste verwendet lineare Indizierung :

Dieser nächste Weg wird nicht schneller sein als der naive Ansatz, aber hier ist eine andere Lösung, um das untere Dreieck mit bsxfun zu berechnen:

Für das obere Dreieck:

Eine andere Lösung für die ganze Matrix mit cumsum für diesen Sonderfall (wo v=1:N ). Dieser ist in der Geschwindigkeit nahe.

Vielleicht können diese ein Ausgangspunkt für etwas Besseres sein.

Die Frage, warum Sie das tun wollen, ist wirklich wichtig.

Im theoretischen Sinne wird der in den anderen Antworten vorgeschlagene Dreiecksansatz Ihnen Operationen ersparen. @ jgmaos Antwort ist besonders interessant beim Reduzieren von Multiplikationen.

Im praktischen Sinne ist die Anzahl der CPU-Operationen nicht mehr die zu minimierende Metrik beim Schreiben von schnellem Code. Speicherbandbreite dominiert, wenn Sie so wenige CPU-Vorgänge haben, so abgestimmt Cache-bewusste Zugriffsmuster sind, wie dies schnell gehen. Der Matrixmultiplikationscode wird extrem effizient implementiert, da es sich um eine so häufige Operation handelt, und jede Implementierung der BLAS-Numerischen Bibliothek, die sich lohnt, verwendet optimierte Zugriffsmuster und auch SIMD-Berechnungen.

Selbst wenn Sie direkt C geschrieben und Ihre Op-Zahl auf das theoretische Minimum reduziert haben, würden Sie wahrscheinlich immer noch nicht die gesamte Matrix multiplizieren. Worauf es ankommt, ist das numerische Primitiv zu finden, das Ihrer Operation am ehesten entspricht.

Alles, was gesagt wurde, es gibt eine BLAS-Operation, die etwas näher kommt als DGEMM (Matrix multiplizieren). Es heißt DSYRK, die Rang-k-Aktualisierung, und es kann für genau A'*A verwendet werden. Die MEX-Funktion, die ich vor langer Zeit geschrieben habe, ist hier. Ich habe es schon lange nicht mehr durcheinander gebracht, aber es funktionierte, als ich es zum ersten Mal schrieb, und lief tatsächlich schneller als eine direkte A'*A .

Dies ist 3 mal schneller als (1: N). * (1: N) vorausgesetzt, ein int32 Ergebnis ist akzeptabel (es ist sogar noch schneller, wenn die Anzahl klein genug ist, um int16 anstelle von int32 zu verwenden ):

Benchmarking:

Beachten Sie, dass aux.'*aux nicht für aux = int32(1:N) verwendet werden kann.

Wie von @ DanielE.Shub gezeigt, muss, wenn das Ergebnis als double -Matrix benötigt wird, eine abschließende Besetzung durchgeführt werden, und in diesem Fall ist die Verstärkung sehr klein: